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MODAL  ACOUSTIC  TRANSMISSION  LOSS  (MOATL): 

A  TRANSMISSION-LOSS  COMPUTER  PROGRAM  USING  A 
NORMAL-MODE  MODEL  OF  THE  ACOUSTIC  FIELD 
IN  THE  OCEAN 


INTRODUCTION 

Investigations  of  acoustic  transmission  in  shallow  water  typically  consider  propagation  to  ranges  of 
100  to  1000  times  the  thickness  of  the  water  column.  The  ocean  may  therefore  be  considered  to  be  a 
thin  film  through  which  the  signal  is  propagated.  In  many  ocean  areas  the  thickness  of  this  film  will 
show  considerable  variation  over  propagation  ranges  of  interest.  In  addition,  the  acoustical  properties 
of  the  water  and  the  ocean  bottom  can  depend  upon  range.  However,  these  range  dependences  are 
usually  slow  and  the  acoustical  properties  are  uniform  over  range  intervals  many  times  the  acoustic 
wavelength.  In  contrast  to  the  slow  range-dependence,  the  acoustical  properties  of  the  ocean  bottom 
and  the  water  column  frequently  show  rapid  depth-dependence.  Appreciable  changes  in  these  media 
occur  over  vertical  distances  comparable  to  the  acoustic  wavelength.  In  addition  to  this  spatial  variabil¬ 
ity,  the  properties  of  the  water  column  are  time-dependent  and  can  change  considerably  over  the  period 
of  a  day  due  both  to  diurnal  heating  and  cooling  and  to  tidal  flow. 

The  relative  shallowness  of  the  water  column  and  the  strong  depth-dependence  of  its  acoustical 
properties  make  normal-mode  representations  of  the  acoustic  field  more  useful  and  reliable  than  the 
ray-tracing  methods  frequently  used  in  transmission  loss  calculations  for  the  deep  ocean.  NRL  has 
developed  a  computer  program  that  calculates  transmission  loss  by  using  a  normal-mode  model.  The 
subroutines  which  perform  the  normal-mode  calculations  have  been  described  previously  [1],  The 
transmission-loss  calculation  using  the  normal-mode  parameters  is  the  subject  of  the  present  report. 

This  transmission-loss  model  may  be  used  with  arbitrary  depth-dependence  of  the  sound  speed  in 
the  water  column  and  in  the  sediment  layer.  Provision  is  made,  via  the  adiabatic  approximation,  for 
calculating  loss  in  an  environment  changing  slowly  with  range.  The  third  source  of  variability  men¬ 
tioned  above,  temporal  change  in  the  water  column,  is  not  considered. 

In  the  following  section  the  normal-mode  model  of  the  acoustic  field  is  described  briefly.  Details 
of  this  model  and  the  associated  FORTRAN  programs  are  found  in  Ref.  1.  Recent  revisions  of  these 
programs  are  described  in  Appendix  A.  The  transmission-loss  model  for  the  coherent  and  incoherent 
modal  field  sum  for  the  perfectly  stratified  (range-independent)  ocean  is  then  presented.  Modifications 
made  to  the  calculation  to  incorporate  an  environment  changing  slowly  with  range  follow. 

THEORY 

Normal-Mode  Model  for  a  Perfectly  Stratified  Medium 

The  model  geometry  is  shown  in  Fig.  I.  A  fluid  layer  of  thickness  H\  and  uniform  density  p\  is 
bounded  above  by  a  pressure-release  surface  and  below  by  a  second  fluid  layer,  which  has  thickness  H: 
and  uniform  density  p2.  These  layers  will  be  referred  to  as  the  water  layer  and  the  sediment  layer.  The 
(arbitrary)  sound  speed  profiles  in  the  water  and  sediment  layers  are  C|  (z)  and  c2  (z),  respectively. 
Beneath  the  sediment  layer  is  a  homogeneous  semi-infinite  basement  of  uniform  density  p2  and 
compressional  sound  speed  c}c.  The  basement  may  be  modeled  as  a  fluid  or  as  a  shear-supporting 
solid.  In  the  latter  case,  the  shear  sound  speed  is  c3s. 
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Fig.  1  —  Physical  model.  An  infinite  half-space  consisting  of  two  fluid  layers 
bounded  above  by  air  and  having  respective  depths  H,  and  H2,  densities  and 
p2,  and  a  third,  semi-infinite  layer  of  density  p3,  compressional  velocity  c3c, 
and  shear  velocity  c3s  (if  it  is  a  solid).  At  the  source  point  the  2-axis  of  a 
cylindrical  coordinate  system  is  established  perpendicular  to  the  pressure-release 
surface  (the  r-axis)  with  increasing  z  downward.  The  sound  speed  profile  in  the 
water  and  sediment  layers,  c3(r)  and  c2(z)  respectively,  is  a  function  of  depth. 

A  cylindrical  coordinate  system  is  defined  so  the  pressure-release  surface  lies  in  the  ( r ,  0)  plane, 
and  the  z-axis  is  perpendicular  to  the  surface  with  z  increasing  downward.  A  harmonic  point  source  of 
unit  strength  and  angular  frequency  w  lies  on  the  z-axis  at  depth  z0.  The  velocity  potential  4>  at  any 
field  point  (r,  0,  z)  satisfies  the  wave  equation: 

V2«t+  4>  —  —  —  8(r)8(0)8(z  —  z0).  (1) 

c(z)  f  r 

The  model  geometry  possesses  cylindrical  symmetry.  The  boundary  conditions  at  the  media  inter¬ 
faces,  the  water  depth,  and  c(z)  do  not  depend  on  r,  so  we  may  separate  Eq.  (1)  into  two  ordinary 
differential  equations.  The  resulting  solution  is: 

*(r,  0-JJ7-  p(Co>  I  ««(Co)  «.({)  » o'”  (k„r). 


where  N  is  the  number  of  discrete  normal  modes  allowed  and  where  we  have  introduced  the  dimen¬ 
sionless  depth  coordinate  {  -  z/H\.  The  eigenfunctions  u„(0  satisify  the  eigenvalue  equation: 
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and  in  the  case  of  a  fluid  basement  are  subject  to  the  normalization  condition 

f0  p(Ou?(Od£=  1, 


(3) 


where,  depending  on  the  value  of  £,  p  —  plt  p2,  or  p3  and  c  =  C\,  c2,  or  c3.  A  similar  (but  more  com¬ 
plicated)  condition  applies  to  the  solid-basement  model  [1J.  These  normalized  eigenfunctions  are  the 
acoustic  normal  modes  of  the  given  environment.  At  sufficiently  long  range  the  Hankel  function 
//o")  ik„  r )  may  be  replaced  by  its  asymptotic  form.  Thus 

1  )l/2  ^  M„(fo )u»(0  Hk„r-iol- ir/4) 


where  the  time  dependence  e~'ul  has  been  inserted. 


.1/2 


e 


(4) 


Each  of  the  terms  in  the  sum  in  Eq.  (4)  corresponds  to  the  contribution  of  a  single  normal  mode 
of  propagation.  Each  of  these  modal  contributions  is  propagated  independently  of  the  others.  Attenua¬ 
tion  of  the  signal  field  is  introduced  by  allowing  the  wave  number  (eigenvalue  of  Eq.  (2))  of  each  mode 
to  become  complex: 

k„~  kn  +  /S„. 

The  attenuation  coefficient,  8„,  assumes  the  form: 

8  „  =  «2Vn2)  +  ejc-Tn3"’  +  fi$y(n3s)  +  $<).«  +  S|.„  +  «„•  <5) 


Here  e2  is  the  plane-wave  attenuation  coefficient  (imaginary  part  of  the  wavenumber)  in  a  hypothetical 
infinite  medium  consisting  of  the  material  in  the  sediment  layer.  The  quantities  e3(.  and  e3s  represent 
the  compressional  and  shear  plane-wave  attenuation  coefficients  of  the  basement.  The  quantities  y(„2\ 


y(n3s)  measure  the  degree  to  which  the  nth  mode  interacts  with  the  sediment  and  the  basement 


compressional  and  shear  wave  mechanisms.  If  the  basement  is  a  fluid,  the  term  e3sy{„is)  is  absent.  Of 
the  remaining  terms,  S0  „  and  5,  „  represent  attenuation  of  the  modal  field  due  to  interaction  of  the 
mode  with  statistically  rough  boundaries  at  the  pressure-release  boundary  and  the  water-sediment 
boundary,  respectively.  The  rough-boundary  interaction  is  discussed  in  Refs.  2  and  3.  The  term  an 
represents  the  attenuation  due  to  absorption  by  the  water  (see  Appendix  A).  The  inclusion  of  the 
attenuation,  Eq.  (5),  due  to  rough  boundaries  and  water  and  bottom  absorption  in  Eq.  (4)  gives: 

l1/2  »  uMo)un(0 


<t»(r,  t)  —  />(< 0) 


1 

8ir//|  r  I 


i(kH 


k'J 2 


uf— tt/4)  —5 

e 


(6) 


The  (real)  instantaneous  pressure  pit)  due  to  a  signal  source  of  rms  source  pressure  level  S,  referred  to 
unit  distance  from  the  source,  is: 


pit)  -  S(4ir)!/2p(Co)  £  T7~  THT  cos iknr  -  cot  -  it/ 4)  e  ^ 


iHxknrV' 

Details  of  the  results  presented  in  this  section  are  given  in  Ref.  4. 


(7) 


Transmission  Loss  for  a  Perfectly  Stratified  Medium 


To  obtain  transmission  loss  we  consider  the  rms  pressure  averaged  over  a  time  T  »  — : 

OJ 

1/2 


<p2it)> 


1/2 


S(4w) 


*n)'12  1  rT 

H{  Pl|  T  Jo 


^  U„iCo)u„iO  „  -#„r| 

S  — 7, — 77/2 —  cos(kn  r  -  cot  -  tt/4)  e 

n-i  \k„r ) 


dt 
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This  expression,  in  which  the  summation  includes  the  phases  of  the  individual  modal  pressure  contri¬ 
butions,  is  called  the  coherent  sum.  The  coherent  transmission  loss  obtained  from  this  expression, 
expressed  in  decibels,  is: 


L,oh  =  -10  log 


(2np\) 


Hi 


y  e~S"'cosUc  r) 

h  U„r)I/2  " 


"  uMo)u„  (Q 
h  U„r),/2  ' 


sin  (kn  r)  |j. 


(8) 


When  N  is  large,  loss  calculated  from  this  expression  usually  exhibits  rapid  oscillations  of  order  10  to 
20  dB  as  range  changes  (see  Fig.  2).  Transmission  loss  measurements  employing  CW  acoustic  signals 
show  similar  oscillations  (see  Fig.  3).  These  oscillations  are  caused  by  phase  interference  effects  among 
the  normal  modes  in  which  the  signals  are  propagated.  Details  of  the  interference'  pattern  are 
extremely  sensitive  to  the  values  of  k„.  The  values  of  k„  are,  in  turn,  sensitive  to  the  sound-speed 
structure  of  the  ocean  bottom.  In  most  cases  of  practical  interest  when  there  are  more  than  a  few 
modes  the  sound-speed  structure  of  the  ocean  bottom  is  not  known  with  sufficient  accuracy  to  permit 
detailed  agreement  between  calculated  and  measured  interference  patterns.  Comparison  of  calculated 
and  measured  results  is  aided  if  the  rapidly  varying  interference  pattern  is  removed,  leaving  only  a 
smooth  curve.  In  treating  experimental  data  this  is  accomplished  by  smoothing  CW  loss  measurements 
over  a  range  interval  or  by  using  broadband  signals  and  processing  techniques.  The  interference  pattern 
is  removed  from  the  model  calculations  by  performing  an  incoherent  mode  summation.  That  is,  the 
energy  contributions  of  individual  modes  rather  than  the  phased  pressures  are  added.  The  resulting 
expression  for  the  incoherent  loss  is 


Linc  =  —  101og10 


(2ttpi2) 

~hT 


N 

I 

n—  1 


"n  ^O)  UJO  -<vl2| 
(*„r),/2  J 


(9) 


Treatment  of  Nearly  Stratified  Media 

In  most  shallow  water  areas  of  interest  the  assumption  that  the  geometry  and  acoustical  properties 
of  the  medium  do  not  depend  upon  range  is  not  valid,  even  over  relatively  short  (—10  km)  propaga¬ 
tion  paths.  When  a  range-dependent  medium  is  introduced,  the  acoustic  wave  equation  (Eq.  (1))  can¬ 
not  be  treated  by  the  separation-of-variables  technique  used  above.  Since  solutions  of  this  generalized 
problem  do  not  exist,  it  is  necessary  to  employ  approximation  techniques.  The  approximation  used 
here  is  that  the  range-dependence  of  the  environment  is  sufficiently  slow  that  the  wave  equation  is 
"locally  separable."  By  this  we  mean  that  any  property  of  a  given  normal  mode,  say  the  eigenvalue  k„  or 
the  attenuation  coefficient  8„,  in  the  vicinity  of  some  point  in  the  range-dependent  medium  is  the  same 
as  it  would  be  in  a  hypothetical  range-independent  medium  with  an  environment  the  same  as  at  the 
point  of  interest.  In  other  words,  the  normal  modes  of  propagation  adapt  to  the  local  environment  and 
the  local  properties  can  be  calculated  from  the  range-independent  model. 

An  additional  approximation,  that  the  range-dependent  environment  does  not  transfer  energy 
from  one  mode  to  another,  is  made.  In  this  approximation  [5,6],  called  the  adiabatic  approximation  or 
the  conservation  of  mode  index,  energy  originally  propagated  in  a  particular  normal  mode  remains  in 
that  mode  until  it  is  removed  by  absorption. 

The  modifications  [7]  to  Eqs.  (8)  and  (9)  necessary  to  employ  these  two  approximations  are: 

(Jo)  U„  (£)  —  u„  (Co)  u'„  (c) 

and 


W,  -  Vtf,  H\ , 


4 


Range  (km) 
(c) 


Fig.  2  —  Calculated  transmission  loss  These  three  graphs  are  examples  or  the  plotted 
output  generated  by  PROGRAM  MOATL.  All  result  from  the  same  physical  environ¬ 
ment  (given  as  'Test  Case  Number  1"  in  (he  "OUTPUT"  section  of  this  report),  but  each 
corresponds  to  a  different  value  of  the  receiver  depth:  (a)  70  m,  (b)  200  m,  and  (c)  400 
m  The  value  of  l.coh  (sec  Eq  (8))  is  plotted  as  a  continuous  line  and  exhibits  the  oscil¬ 
lations  due  to  modal  interference  The  value  of  L,„r  at  each  range  (see  Eq.  (9))  is  plotted 
as  a  circle  The  circles  overlap  at  most  of  the  ranges  in  these  illustrations 


MILLER  AND  WOLF 


RRNGE  [KM] 

Fig.  3  —  Measured  CW  transmission  loss.  Oscillations  in  the  transmission  loss  are  shown 
here  for  some  typical  measurements  employing  a  towed  CW  source.  The  places  where 
data  are  missing  correspond  to  intervals  during  which  the  source  was  turned  off. 


where  unprimed  quantities  u„  and  H\  apply  to  the  source  location  and  primed  quantities  w„  and  H\ 
apply  to  the  field  point  at  range  r  ; 


is  the  cumulative  phase;  and 


k„r  =  fQ  k„(r)  dr 

^  Jo  Bn  ( r )  dr 


is  the  average  attenuation  coefficient.  Eqs.  (8)  and  (9)  then  become; 

(27rp  ,2) 
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The  phase  of  the  signal  is  obtained  as  the  arctangent  of  the  ratio: 


N 

u„(io)u^0 

-A_r  . 

I 

n—  1 

-  » 

<Pn 

e  sin0„ 

N 

un(U)un(0 

-hnr  , 

I 

n-1 

<Pn 

e  cosd>„ 

(10) 


(11) 


ANALYSIS  OF  PROGRAM  AND  TECHNIQUES 
General  Remarks 

For  each  point  source  of  harmonic  frequency  /  (/=«/2ir )  PROGRAM  MOATL  calculates  the 
transmission  loss  as  a  function  of  range  r  and  depth  z.  The  model  described  above  located  the  source  at 
the  range  origin  and  considered  the  field  point  associated  with  the  receiver  to  be  a  variable.  In  the  pro¬ 
gram  code,  however,  source-receiver  reciprocity  is  employed  to  locate  the  receiver  at  the  range  origin 
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and  depth  z0.  The  source  is  located  at  the  point  (r,z).  This  change  was  made  so  that  the  program  con¬ 
forms  to  the  common  experimental  situation  in  which  a  source  is  lowed  along  a  radial  track  extending 
from  a  fixed  receiver.  In  model  calculations  for  range-independent  environments  the  location  of  the 
receiver  at  the  origin  is  largely  for  convenience  of  notation.  In  range-dependent  environments,  how¬ 
ever,  locating  the  fixed  transducer  at  the  origin  makes  the  calculation  more  economical,  since  <f>„  and 
A„  may  be  evaluated  over  the  track  without  repetitive  integrations  over  range.  The  program  described 
in  this  report  may  be  applied  to  situations  in  which  the  source  is  fixed  and  the  receivers  are  moving,  if 
the  user  exchanges  the  variables  associated  with  source  and  receiver. 

In  the  following  description  and  in  programming  Eqs.  (10)  and  (11),  «„(£ 0)  is  taken  to  be  the  «th 
modal  eigenfunction  calculated  by  using  the  environment  found  at  the  origin  and  evaluated  at  the  (nor¬ 
malized)  depth  £o  associated  with  the  receiver.  The  quantity  «„(£)  is  calculated  for  the  environment 
associated  with  range  r  and  is  evaluated  at  the  source  depth  £.  Note  that  if  the  water  depth  at  the 
source  H{  depends  upon  range,  the  normalized  source  depth  £(/•)  =  z/ H[  ( r )  is  also  range-dependent. 

The  needed  values  u„  are  obtained  from  the  normal-mode  representation  of  the  sound  field  in  an 
ocean  consisting  of  a  three-layer  half-space,  as  discussed  previously.  To  this  end,  the  main  program 
(MOATL)  calls  one  of  two  sets  of  subroutines;  the  set  FLUID,  HALFF,  and  ITRTF  assume  the  semi¬ 
infinite  basement  to  be  a  fluid,  whereas  the  set  SOLID,  HALFS,  and  ITRTS  assume  a  basement  capable 
of  supporting  shear.  These  subroutines  are  streamlined  versions  of  the  programs  described  in  Ref.  1. 
Some  improvements  have  been  made  since  Ref.  1  was  published,  the  major  changes  are  discussed  in 
Appendix  A.  Variable  namelists  are  given  in  Appendix  B. 

When  two  distinct  bottom  layers  with  different  properties  are  not  required,  a  two-layer  model  may 
be  realized.  This  is  accomplished  by  giving  the  sediment  layer  (a)  a  small  thickness  and  (b)  physical 
properties  identical  to  those  either  of  the  water  layer  (SOLID  or  FLUID)  or  of  the  basement  (FLUID 
only).  Realization  of  the  single-medium  bottom  is  illustrated  by  the  first  test  case  in  the  "OUTPUT" 
section. 

Provision  is  made  for  the  receiver  to  be  located  in  the  water  or  sediment  layer;  the  source  must 
be  in  the  water  layer. 

The  range  points  at  which  loss  is  to  be  calculated  or  the  calculation  ranges,  as  we  shall  sometimes 
call  them,  are  assumed  to  be  equally  spaced.  The  user  supplies  the  number  of  points  and  the  max¬ 
imum  range  at  which  calculated  transmission  loss  is  desired.  The  program  obtains  the  spacing  between 
points  by  division.  Nevertheless,  it  is  a  simple  matter  to  obtain  results  ?.i  unequally  spaced  range 
points.  There  are  two  places  in  the  main  program  where  the  FORTRAN  code  must  be  altered  slightly. 
These  places  are  marked  with  COMMENT  statements  which  give  examples  of  how  to  make  the 
required  modifications  (see  the  listing  given  in  Appendix  C).  One  word  of  caution;  difficulties  may  be 
encountered  if  NMFREQ  is  given  a  value  other  than  one;  unequally  spaced-calculation-range  cases 
should  be  run  one-at-a-time. 

The  user  must  always  supply  an  environmental  data  set  (including  water  depth,  sound-speed 
profile,  bottom  properties,  and  other  relevant  parameters)  at  zero  range.  If  a  range-independent  model 
is  to  be  used,  no  further  environmental  input  is  required.  For  a  range-dependent  model  calculation, 
additional  environmental  data  sets  are  required.  Each  set  may  differ  from  the  others  by  any  arbitrary 
combination  of  the  environmental  parameters,  subject  to  the  constraint  that  the  fluid-basement  and 
solid-basement  models  may  not  both  be  used  along  the  track.  Such  data  sets  may  be  supplied  at  any 
number  of  user-selected  ranges  (not  necessarily  equally  spaced  and  not  necessarily  coincident  with  any 
of  the  ranges  at  which  loss  is  to  be  calculated) .  Guidelines  for  the  selection  of  spacing  of  environmen¬ 
tal  data  sets  are  given  in  the  "INPUT  DATA"  section  below.  The  normal-mode  parameters  are  calcu¬ 
lated  at  each  of  these  ranges.  The  modal  properties  are  obtained  at  the  intermediate  range  points 
(where  the  loss  is  to  be  calculated)  by  performing  a  linear  interpolation.  If  loss  calculations  are  desired 
out  past  the  last  data  set,  a  linear  extrapolation  is  performed. 
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With  regard  to  this  interpolation,  there  are  four  range  points  of  interest  at  any  given  moment  dur¬ 
ing  the  execution  of  the  program  (for  a  range-dependent  model).  The  first  is  the  range  origin.  Since 
an  environmental  profile  is  always  supplied  here,  no  range-interpolation  is  necessary  to  obtain  the 
modal  parameters  at  the  receiver  position.  The  second  range  of  interest  is  the  calculation  range  /•, 
which  together  with  the  source  depth  locates  the  source  position.  Since  r  is  not  usually  coincident  with 
the  range  of  one  of  the  user-supplied  input  environments,  a  linear  interpolation  between  the  two  input 
ranges  which  bound  r  is  required  in  order  to  obtain  the  modal  parameters  at  the  source  position.  The 
ranges  at  which  environmental  data  are  supplied  are  the  third  and  fourth  ranges  of  interest.  The 
smaller  input  range  or  any  group  of  variables  at  that  range  will  be  designated  hereafter  by  SR;  the  larger 
input  range  by  LR.  Details  of  the  interpolation  are  presented  below  in  Part  7  of  the  "Step-by-Step 
Analysis"  section. 

The  dimensioning  of  many  of  the  arrays  implies  an  application  for  which  the  number  of  normal 
modes  is  less  than  or  equal  to  150.  If  more  modes  are  expected,  redimensioning  and  some  minor 
FORTRAN  code  modifications  are  required;  these  are  not  discussed  here.  Knowledge  of  the  expected 
number  of  modes  is  also  required  for  determination  of  the  input  variables  Lll  and  LI2.  For  this  pur¬ 
pose,  we  give  the  following  "rule  of  thumb"  guide; 

_/T~—  J+  1/2,  (12) 

where  N  is  the  total  number  of  modes,  /  is  the  frequency,  H  is  the  sum  of  the  thicknesses  of  the  water 
and  sediment  layers,  cmin  is  the  minimum  sound  speed  found  in  the  water  and  sediment  layers,  and  c} 
is  c3t.  (FLUID  case)  or  c3s  (SOLID  case).  For  an  isovelocity  profile  and  slight  redefinition  of  variables, 
the  expression  becomes  exact  [81. 

Before  using  the  program,  one  must  first  inspect  the  two  PARAMETER  statements  at  the  begin¬ 
ning  of  the  main  program  and  change  them  if  necessary.  The  variables  REC  and  SOC  should  be 
assigned  values  equal  to  or  greater  than  the  number  of  receiver  and  source  depth  points,  respectively, 
at  which  calculated  results  are  desired.  The  variable  REC5  should  be  equal  to  or  greater  than  REC  and 
should  be  an  integral  multiple  of  five.  This  parameter  is  used  for  dimensioning  the  output  transmission 
loss  arrays,  but  the  purpose  of  introducing  it  in  place  of  REC  is  solely  to  make  the  printout  more 
aesthetic;  if  resulting  program  storage  requirements  exceed  the  limitations  of  the  given  computer,  the 
parameter  may  be  eliminated  by  minor  output  adjustments.  The  variable  RNG  should  be  assigned  a 
value  equal  to  or  greater  than  the  number  of  range  points  at  which  calculations  are  desired.  Before 
using  the  program  for  the  first  time  on  a  given  machine,  the  parameters  MGNTD  and  PRCSN  should 
be  assigned  appropriate  values  (see  the  COMMENT  statement  preceding  the  PARAMETER  statement 
in  the  listing).  They  will  not  have  to  be  changed  subsequently. 

The  program  was  written  in  ASC  FORTRAN  for  use  with  the  Texas  Instruments  Advanced 
Scientific  Computer  (ASC)  located  at  NRL.  Wherever  possible,  however,  the  source  code  was  put  into 
standard  form.  Thus  it  should  compile  on  most  FORTRAN  compilers  with  a  minimum  amount  of  pre¬ 
liminary  code-changing. 

The  ASC  has  a  single  precision  floating  point  word  (32  bits)  consisting  of  1  bit  for  the  sign,  7  bits 
for  the  exponent,  and  24  bits  for  the  fraction  (precise  to  approximately  7  decimal  digits).  Some  of  the 
program  variables  are  in  DOUBLE  PRECISION.  A  double  precision  word  (64  bits)  consists  of  1  bit  for 
the  sign,  7  bits  for  the  exponent,  and  56  bits  for  the  fraction  (precise  to  approximately  16  decimal 
digits).  In  general,  any  variable  involved  in  or  affecting  the  calculation  of  an  eigenvalue  or  eigenfunc¬ 
tion  is  in  DOUBLE  PRECISION. 

The  required  storage  allocations  for  the  main  program,  subroutines,  and  COMMON  blocks  are 
given  in  Table  1. 
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Table  1  —  Storage  Requirements 


Routines 

Number  of  Words 
(in  hexadecimal  base) 

MOATL 

4C59 

FLUID 

195B 

1TRTF 

4E3 

HALFF 

306 

SOLID 

1B9E 

ITRTS 

581 

HALFS 

335 

Common  Blocks 

TNIH 

2 

TNH 

12C 

TNI 

7 

TH 

1 

TN 

58680 

NIH 

968 

NH 

S 

N1FLU 

2 

NI 

12DB 

IH 

2 

N1SOL 

7 

As  an  aid  to  following  the  flow  of  control  in  the  program  when  reading  the  code,  the  following 
types  of  control  statements  have  been  indented  three  spaces:  (1)  DO  loop;  (2)  GO  TO  statement;  (3) 
transfer-of-control  IF  statement;  and  (4)  calls  to  subroutines. 

The  Naval  Research  Laboratory’s  computer  peripherals  include  an  11-inch  Calcomp  (California 
Computer  Products)  Model  565  Drum  Plotter.  The  on-line  plotter  software  on  NRL's  ASC  currer..l> 
supports  this  plotter.  PROGRAM  MOATL  includes  an  option  which  uses  this  package  to  plot  coherent 
and  incoherent  transmission  loss  as  functions  of  range.  Separate  plots  are  generated  for  each  fre¬ 
quency,  receiver  depth,  and  source  depth. 

Input  to  the  program  is  from  logical  unit  five  (card  reader  by  convention),  and  printed  output  is 
to  logical  unit  six  (line  printer  by  convention). 

Step-by-Step  Analysis 

In  the  present  discussion,  we  follow  PROGRAM  MOATL  step-by-step  from  start  to  finish.  The 
FORTRAN  namelist  (Appendix  B)  should  prove  useful  to  the  reader  at  this  time.  A  synopsis  of  the 
workings  of  the  program  will  be  sketched  and,  whenever  appropriate,  the  numerical  methods  and  pro¬ 
gramming  techniques  will  be  described.  Wherever  the  normal-mode  subroutines  are  mentioned,  refer¬ 
ence  to  Appendix  A  may  prove  useful. 

The  program  has  been  broken  up  into  1 1  parts  for  discussion.  The  program  listing  is  given  in 
Appendix  C.  The  parts  are  defined,  by  control  statement  number  (CSN),  as  follows: 

Parti:  34-58  Part  7:  189-214 

Part  2:  59-66  Part  8:  215-247 

Part  3:  67-120  Part  9:  248-262 

Part  4:  121-141  Part  10:  263-271 

Part  5:  142-146  Part  11:  272-299 

Part  6:  147-188 


9 


MILLER  AND  WOLF 


Preceding  the  executable  code  are  CSNs  1-33,  which  set  up  the  necessary  COMMON  blocks, 
PARAMETERS,  DIMENSIONS,  and  FORMATS. 

Part  1:  Input  Data  and  Initialization 

The  plot  package  is  initialized  and  the  parameters  PMGTD  and  PPRCN  are  defined.  These  two 
parameters  are  used  in  the  normal-mode  subroutines.  The  transmission-loss  input  parameters  are  next 
read  in  and  printed  out.  The  variable  DR  is  the  calculation-range  increment,  defined  by  dividing  the 
maximum  range  by  the  number  of  calculation-range  points.  The  array  RAN(l)  contains  integral  multi¬ 
ples  of  DR,  which  are  the  ranges  at  which  calculations  are  to  be  performed.  If  plotting  is  desired, 
parameters  are  now  defined  for  this  purpose. 

Part  2:  Frequency  and  Attenuation 

The  frequency  loop  (DO  340)  marks  the  beginning  of  an  actual  calculation  of  transmission  loss. 
A  value  of  NMFREQ  greater  than  one  may  be  used,  not  only  to  obtain  results  for  more  than  one  fre¬ 
quency,  but  for  more  than  one  run  of  the  program  for  any  reason  (e.g.  different  sediment  thickness, 
different  profile,  etc.).  For  each  case,  the  TITLE  and  frequency  F  are  read  in  and  printed  out.  The 
equation  for  EP4  converts  the  plane-wave  absorption  coefficient  of  the  sediment,  e2  (EP1  in  the  FOR¬ 
TRAN),  from  units  of  dB/Hz-m  into  units  of  nepers/m,  in  which  form  it  is  subsequently  used.  Similar 
equations  convert  the  basement  compressional  plane-wave  absorption  coefficient,  ej,  (EP2  in  the  FOR¬ 
TRAN),  and  the  basement  shear  plane-wave  absorption  coefficient,  eis  (EP3  in  the  FORTRAN).  They 
become  EP5  and  EP6,  respectively. 

Part  3:  Input  Data,  Initialization,  Receiver  Parameters 

The  environmental  input  parameters  (at  zero  range)  needed  for  the  normal-mode  calculations  at 
the  site  of  the  receiver  are  first  read  in  and  printed  out.  The  appropriate  normal-mode  subroutines  are 
next  called  to  perform  modal  calculations.  Prior  to  the  call  to  FLUID  or  SOLID,  NMODE  is  set  to 
10000.  This  is  done  for  the  following  reason.  In  a  range-dependent  calculation,  one  of  the  input 
environments  may  support  more  modes  than  a  previous  environment,  i.e.,  one  at  closer  range.  How¬ 
ever,  the  program  implements  conservation  of  mode  index  by  excluding  the  higher  order  modes  which 
are  not  present  at  the  previous  environment.  For  example,  if  only  five  modes  exist  at  an  input  range  of 
10  km,  then  at  each  of  the  calculation  ranges  beyond  10  km,  only  the  five  lowest  order  modes  will  be 
used  for  a  calculation  of  transmission  loss.  Additional  modes  allowed  at  ranges  greater  than  10  km  are 
assumed  to  be  cut  off  at  10-km  range.  Each  time  FLUID  or  SOLID  is  called  at  a  new  input  range, 
NMODE  is  redefined  to  be  the  smaller  of  (a)  the  previous  value  of  NMODE  or  (b)  the  maximum 
number  of  modes  existing  for  the  given  environment.  Since  this  test  is  performed  even  for  the  first  call 
to  FLUID  or  SOLID,  NMODE  must  have  been  defined  prior  to  the  first  subroutine  call.  Since 
"redefinition"  is  actually  to  be  definition  by  criterion  (b),  NMODE  must  be  preset  to  a  large  number. 

Prior  to  the  first  call  to  FLUID  or  SOLID,  KA  is  set  to  zero.  The  variable  KA  is  a  flag  which 
when  zero  causes  the  eigenfunctions  to  be  stored  in  UNRMKIM,!)  and  when  one  causes  the  eigen¬ 
functions  to  be  stored  in  UNRM2(IM,1).  Mode  order  is  designated  by  the  variable  IM,  depth  index  by 
the  variable  I. 

Many  of  the  variables  defined  in  Part  3  have  names  ending  with  the  numeral  1  or  2,  for  example 
RANGE1  and  RANGE2.  The  reason  for  this  (the  same  as  for  UNRM1(1M,I)  and  UNRM2(IM,I)) 
rests  in  the  numerical  technique  employed  to  calculate  the  transmission  loss  for  a  range-dependent 
environment.  Any  given  source  range  at  which  loss  calculations  are  desired  will  fall  between  two 
ranges  at  which  environmental  input  data  have  been  supplied.  Environmental  and  modal  parameters 
required  at  the  source  position  are  approximated  by  a  linear  interpolation  which  uses  the  given  input  at 
each  of  the  bounding  ranges.  (This  procedure  is  described  later.)  The  parameters  at  the  smaller  (i.e.. 
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closer  to  the  receiver)  input  range  (hereafter  designated  SR)  are  stored  in  the  variables  whose  names 
contain  the  trailing  numeral  1.  The  parameters  at  the  larger  range  (LR)  use  the  trailing  numeral  2. 
Parameters  used  for  interpolation  are.  HU,  H12  (water  layer  thickness);  RANGE1,  RANGE2  (range 
at  which  environmental  input  is  supplied);  CT1,  CT2  (sound  speed  at  the  surface  of  the  water  layer); 
CB1,  CB2  (sound  speed  at  the  bottom  of  the  water  layer);  N11,N12(N1  =  LI1  +  1,  where  Lll  is  the 
number  of  incremental  intervals  into  which  the  water  layer  is  broken);  E1GVLK1M),  EIGVL2HM) 
(the  eigenvalue  /c„);  R1H1M),  R12(1M)  (sediment  attenuation  ratio  y{„2))\  R21GM),  R220M)  (base¬ 
ment  compressional  attenuation  ratio  -y(S3^',);  R310M),  R32(1M)  (basement  shear  attenuation  ratio 
RAl(lM),  RA2(IM)  (water  absorption  a„  —  see  Appendix  A);  RTK1M),  RT20M)  (r0,n  — 
see  below);  RBl(IM),  RB2GM)  (T  i  „  —  see  below). 

The  quantities  described  above  are  read  in  and/or  calculated  in  Part  3  for  the  receiver  (zero 
range);  they  therefore  constitute  the  first  SR  group.  They  are  stored  in  the  variables  designated  by  the 
trailing  numeral  1.  They  are  also  initially  stored  in  the  LR  group,  for  a  reason  which  is  explained  in  the 
discussion  of  Part  6.  The  terms  S0  „  and  Sln  appearing  in  Eq.  (5)  may  be  rewritten  as: 

•So, n  =  (1  ~  I  T 0 ,n 

S\.„  =  (i  -  !/?,.„!) r,.„  •  (13) 

The  variable  R o,„  is  the  plane-wave  reflection  coefficient  at  the  air/ water  interface,  and  Rt  „  is  the 
plane-wave  reflection  coefficient  at  the  interface  between  the  water  and  sediment  layers.  T0  n  and  T|  „ 
are  the  respective  scattering  ratios,  which  are  calculated  in  the  normal-mode  subroutines.  (See  also 
Ref.  1.)  At  the  receiver  site  (RANGE1  =  0)  these  quantities  are  stored  in  the  SR  variables  RTl(IM) 
and  RBl(lM),  respectively. 

The  quantities  H10,  RERHOl,  and  RERH02  are  defined  as  the  water  depth,  water  density,  and 
sediment  density,  respectively,  at  the  receiver.  They  are  defined  because  H 1 1 ,  RHOl,  and  RH02  will 
take  on  new  values  when  subsequent  input  environments  are  read  in;  however,  the  values  of  these 
quantities  at  the  receiver  will  be  needed  in  the  final  transmission-loss  calculation. 

The  arrays  SE(IM),  S1«M),  S2UM),  S3UM),  SA(1M),  ST(1M),  and  SB(IM)  will  be  described  in 
Part  8.  They  are  now  initialized  to  zero. 

The  final  calculations  of  Part  3  are  to  obtain  the  values  of  u„(£o)>  w,1'cfi  appear  in  Eqs.  (10)  and 
(11).  The  outer  DO  110  loop  varies  the  mode-order  index  n  (programmed  as  IM)  from  1  to  NMODE. 
The  inner  DO  100  loop  varies  the  receiver  identification  index  (programmed  as  Jl)  from  1  to  NDRE 
(the  total  number  of  receiver  depths  supplied).  Each  of  the  receiver  depths  corresponds  to  a  different 
value  of  £0  —  zo/^i-  The  quantity  un((0)  is  programmed  as  RE  (Jl,  IM). 

The  eigenfunction  for  a  given  mode  is  calculated  in  the  normal-mode  subroutines;  values  of  the 
function  are  defined  at  each  of  the  N1  +  N2  —  (Lll  +  1)  +  (LI2  +  1)  incremental  depths  (see 
Appendix  A  and  Ref.  1).  The  receiver  depth,  however,  will  generally  lie  between  two  of  these  incre¬ 
mental  depths.  The  program  performs  a  linear  interpolation,  as  follows,  to  obtain  the  value  of  un({0). 
Assume  for  convenience  that  the  receiver  is  in  the  water  layer;  the  calculations  for  a  receiver  in  the 
sediment  are  similar.  The  program  first  defines  A1  to  contain  the  number,  plus  fraction,  of  incremental 
layers  (numbered  downward  from  the  air/water  surface)  which  corresponds  to  the  receiver  depth.  For 
example,  if  the  receiver  is  exactly  in  the  middle  of  the  third  incremental  layer,  Al  =  2.5.  The  term  IA1 
contains  the  (truncated)  integer  value  of  Al;  following  the  above  example,  1A1  =  2.  Thus  in  general, 
iiAi+i  <  Co  <  Gai+2>  where  Gai+i  is  the  normalized  depth  at  the  top  of  the  (IAl  +  l)ih  incremental 
layer;  Gai+2  is  defined  similarly.  Note  that  G  “  0.  In  the  above  example,  £„  is  bounded  by  the  depths 
at  the  tops  of  the  third  and  fourth  incremental  layers.  (The  general  procedure  is  illustrated  in  Fig.  4.) 
Standard  linear  interpolation  yields  the  value  for  i/„($0): 

W,((o)  “  U„(Ga1  +  |)  +  Mh„({|AI+2)  —  U„(Ga)  +  |)1-  (14) 
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Fig.  4  —  Receiver  eigenfunction.  The  first  mode  (fluid-basement  model)  at 
zero  range  is  illustrated,  along  with  the  values  of  the  eigenfunction  to  be  used 
in  the  interpolation  of  u„(( 0). 

The  quantity  A  (programmed  as  DLTAl)  is  that  fraction  of  a  layer  increment  by  which  f 0  exceeds 
£ui+i-  For  a  receiver  in  the  water  layer,  we  have: 

A  =  (LID  *  (Co)  -  1FIX[(LI1)  *  (Co))  -  -Co~-Cj71-1  ■  (15) 

tlAI  +  2  -  MAl  +  t 

For  a  receiver  in  the  sediment,  the  calculation  of  A  is  performed  similarly. 

For  a  receiver  depth  equal  to  the  depth  of  one  of  the  incremental  layer  boundaries,  A  will  take  on 
the  value  zero  or  one,  and  the  interpolation  is  actually  a  "do  nothing"  procedure. 

Part  4:  Range-Independent  Parameters 

If  the  calculation  is  to  be  based  on  a  range-independent  model,  then  the  transmission- loss  param¬ 
eters  and  the  calculated  eigenfunctions  at  any  given  range  will  be  the  same  as  those  already  calculated  at 
zero  range  for  the  receiver.  Later  calculations  of  the  program  use  the  terms  HS,  CTS,  and  CBS  for  the 
water-layer  thickness  and  sound  speeds  at  the  surface  and  bottom  of  the  water  layer,  respectively.  For 
a  range-independent  calculation,  these  parameters  and  the  others  defined  in  Part  4  will  not  depend  on 
the  range  r  as  the  calculation  point  moves  out  in  range.  Thus  they  have  simple  definitions.  The 
definitions  for  a  range-dependent  calculation  are  given  in  Part  7. 

The  procedure  for  obtaining  the  quantities  w„’(£),  which  are  the  eigenfunction  values  at  the 
source,  is  identical  to  the  procedure  described  in  Part  3  for  the  receiver.  As  before,  there  are  two 
loops,  one  for  the  NDSO  source  depths  £  and  one  for  the  n  =  1,  . . .  ,  NMODE  modes.  The  variable 
u„(£)  is  programmed  as  SM(J2,1M),  where  J2  is  the  depth  index  and  IM  the  mode-order  index.  The 
only  difference  in  procedure  between  the  source  and  receiver  calculations  is  that  the  source  must  be  in 
the  water  layer;  i.e.,  it  may  not  be  located  in  the  sediment. 
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The  quantities  WN(IM),  Gl(IM),  G2(!M),  G3(1M),  GAUM),  GT(IM),  and  GB(IM)  have  sim¬ 
ple  definitions,  thus  obvious  meanings,  for  the  range-independent  case  (see  Appendix  B). 

Part  5:  Calculation  Range  Loop 

The  DO  300  loop  uses  I  as  the  index  for  the  NRCALC  calculation-range  points,  which  are  stored 
in  RAN(I)  one-by-one  as  encountered.  If  the  calculations  are  to  be  range-independent,  then  the 
parameters  of  interest  in  Parts  6  and  7  have  already  been  defined  in  Part  4,  and  Part  5  now  transfers 
control  to  Part  8.  This  is  programmed  as.  lF(NRBUF.EQ.I)  GO  TO  209. 

If  the  calculations  are  to  be  range-dependent,  then  two  more  checks  are  made.  Before  calcula¬ 
tions  can  be  performed,  we  require  various  environmental  and  modal  parameters,  which  are  to  be 
obtained  for  a  given  RAN(I)  by  interpolating  these  same  parameters  between  their  known  (or 
subroutine-calculated)  values  at  the  SR  and  the  LR  ranges  (see  Part  3  for  definitions  and  Part  7  for  the 
technique).  Prior  to  the  first  time  through  the  DO  300  loop,  only  the  zero-range  (receiver)  parameters 
have  been  obtained,  and  they  have  been  stored  in  the  SR  group.  The  first  time  through  the  DO  300 
loop  (and  only  the  first  time),  RANGE2  will  be  equal  to  RANGE1,  which  is  equal  to  zero  (see  Part  3). 
If  they  are  equal,  then  Part  5  transfers  control  to  Part  6,  where  the  LR  group  is  established. 

The  second  check  is  made  for  subsequent  loops  over  range.  If,  for  a  given  value  of  I, 
r  =  RAN(I)  lies  between  the  present  values  of  RANGE1  and  RANGE2,  then  an  interpolation  between 
the  present  SR  and  LR  groups  can  and  should  be  made;  Part  5  thus  transfers  control  to  Part  7.  If 
RAN(l),  which  has  just  been  obtained  by  adding  DR  to  the  previous  calculation  range,  is  greater  than 
RANGE2,  then  the  value  of  RANGE2  must  become  the  next  value  of  RANGE1,  and  a  new  RANGE2 
(and  a  new  LR  group)  is  needed  before  the  calculations  can  proceed  (see  Part  6).  Part  5  transfers  con¬ 
trol  to  Part  6.  The  quantity  RANGE2  +  (DR/2)  is  actually  used  for  this  check,  since  for  RANGE2  < 
RAN(I)  <  RANGE2  +  (DR/2),  it  is  more  accurate  to  extrapolate  past  the  present  RANGE2  than  to 
make  the  redefinitions  of  Part  6  and  interpolate  between  subsequent  SR  and  LR  groups.  (See  Fig.  5.) 

Occasionally  (when  the  option  of  unequally  spaced  calculation  range  points  is  used— see  the  "Gen¬ 
eral  Remarks"  section  and  the  COMMENT  statements  in  the  program  listing  following  CS Ns  51  and 
142),  RAN(I)  may  be  larger  than  RMAX,  the  range  of  the  last  input  environment.  In  this  case,  linear 
extrapolation  is  to  be  performed  past  RMAX,  which  is  RANGE2  at  this  point. 

The  formulae  for  extrapolation  are  identical  to  those  for  interpolation  and  are  not  programmed 
separately.  We  speak  below  only  of  interpolation,  but  the  "double  usage"  is  intended. 

Part  ft:  Range-Dependent  Parameters:  Initialization 

As  alluded  to  previously,  this  part  of  the  program  is  entered  only  if  the  present  value  of  RAN(I) 
lies  outside  the  interpolation  interval  (RANGE1,  RANGE2  +  DR/2)  and  thus  the  interval  must  be 
redefined.  To  this  end,  all  of  the  present  LR  group  variables  are  stored  in  the  SR  group,  i.e.,  RANGE2 
becomes  the  new  RANGEl,  and  all  of  the  LR  parameters  become  the  new  SR  parameters.  For  exam¬ 
ple,  we  encounter  FORTRAN  statements  like  CT1=CT2  in  Part  6.  Note  that  the  previous  values  of 
the  SR  variables  are  lost.  (They  will  no  longer  be  needed  for  calculation  of  modal  parameters  at  the 
source  position.)  The  new  RANGE2  and  its  environmental  parameters  are  next  read  in  and  printed  out, 
and  the  new  LR  group  is  established.  Since  the  normal-mode  subroutines  have  been  called  previously 
to  provide  the  modal  parameters  at  zero  range,  the  initial  zero  value  of  K.A  is  now  changed  to  one. 
When  either  FLUID  or  SOLID  is  subsequently  called  to  perform  the  modal  calculations,  the  value 
KA=1  ensures  that  the  eigenfunctions  are  stored  in  the  LR  array  UNRM2(N,K). 

In  Part  3  we  remarked  that  the  receiver  parameters  were  initially  stored  not  only  in  the  SR  group, 
but  also  in  the  LR  group.  The  reason  for  this  now  becomes  apparent,  in  view  of  the  procedure 
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described  above  of  shifting  the  values  of  the  LR  variables  into  the  SR  variables  each  time  a  new  user- 
supplied  environment  is  encountered.  During  the  first  pass  through  the  DO  300  loop  (1  =  1),  the  first 
nonzero  range  environment  is  read  in.  Just  prior  to  storing  it  in  the  appropriate  LR  variables,  however, 
all  of  the  LR  group  is  shifted  into  the  SR  group.  By  predefining  these  two  groups  to  be  identical,  the 
receiver  parameters  are  not  lost  but  are  still  retained  in  the  SR  group.  The  new  LR  group  is  then  esta¬ 
blished.  Following  the  example  of  the  last  paragraph,  we  encounter  such  FORTRAN  statements  as 
CT2—CH1)  in  setting  up  the  new  LR  group.  (The  trailing  numeral  I  in  the  array  name  here  desig¬ 
nates  the  water  layer,  not  the  SR  group.)  Storage  of  the  eigenfunctions  into  the  SR  and  LR  arrays  are 
handled  in  a  different  manner.  The  flag  KA  is  set  to  zero  to  ensure  that  the  receiver  eigenfunctions  are 
stored  in  UNRM1(N,K)  and  KA  is  then  set  to  one  to  ensure  that  the  eigenfunctions  at  the  first 
nonzero  range  are  stored  in  UNRM2(N,K).  (See  above  and  also  Part  3.) 

Part  7:  Range-Dependent  Parameters:  Calculation  of  Source  Paramet's 

Part  7  and  Part  4  perform  the  same  function  and  correspond  to  the  range-dependent  and  range- 
independent  cases,  respectively.  The  quantities  HS,  CTS,  and  CBS  depend  on  r  for  a  range-dependent 
calculation.  For  the  thickness  of  the  water  layer  (the  other  two  quantities  are  defined  similarly),  stan¬ 
dard  linear  interpolation  requires  us  to  put: 

HS  =  Hll  +  A  *  (H12-  Hll), 

where 


.  RAN  (1)  -  RANGE  1  , 

RANGE  2- RANGE  1 

The  quantity  A  is  programmed  as  SCaLE.  For  interpolation  of  the  eigenvalues  and  attenuation  and 
scattering  ratios,  it  is  more  accurate  to  replace  the  numerator  in  Eq.  (16)  by  [RAN(l)  -  (DR/2)] 
-  RANGE1,  which  is  effected  by  the  FORTRAN  statement. 

SCALE  =  SCALE  -  DR/(2.0*  (RANGE2  -  RANGE  D).  (17) 


The  reason  for  this,  illustrated  in  Fig.  6,  is  as  follows. 


Fig  6  —  Approximation  of  Jn  yj2'  ( i  ')dr  The  "piece"  of  ihc  integral 
/ f an )/'  I >  yj21  (r')rfr'  is  approximated  by  I2R  •  y^21  (r"l.  where  r"  is 
equal  io  RAN(I)  -  ( OR/ 2 ) .  The  result  is  shown  wilh  solid  lines,  and  is 
more  accuraic  than  Ihe  result  (shown  with  dashed  lines)  corresponding  io 
r  -  RAN(I) 
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Equations  (10)  and  (11)  for  the  transmission  loss  each  contain  the  term  e  where 

±n  =  l  f'Sn(r’)dr'  (18) 

is  the  average  attenuation  coefficient.  This  integral  may  be  broken  into  a  number  of  separate  integrals 
through  use  of  Eq.  (S).  For  example,  one  term  of  Eq.  (18),  say  will  be 

^n=llJ'y^(r')dr'.  (19) 


At  any  given  range  r  =  RAN(l),  the  integral  is  approximated  by  adding  a  new  "piece,"  representing  the 
integral  from  RAN(I-l)  to  RAN(I),  to  the  stored  value  of  the  (approximated)  integral  from  zero  to 
RAN(l-l).  (See  Part  8.)  Reference  to  Fig.  6  demonstrates  that  the  appropriate  value  of  y[n2)  to  be 
multiplied  by  DR  (in  order  to  approximate  /ran/i-d  ?i2>  (r')  dr'  by  a  rectangle)  is  the  value  of  the 
function  at  the  midpoint  of  the  range-increment  interval,  i.e.,  at  RAN(l)  -  (DR/2).  In  the  FORTRAN 
code,  this  value  of  y(„2)  is  programmed  as  Gl(IM).  Thus  if  we  write 

G  KIM )  *  Rll(lM)  +  SCALE  *  (R12(1M)  -  Rll(IM)), 


SCALE  must  be  given  by  its  redefined  value,  Eq.  (17),  rather  than  its  original  value,  Eq.  (16). 

Other  variables  handled  in  a  manner  identical  to  that  of  Gl(lM)  are.  WN(IM)  (the  eigenvalue 
kn ),  G2(1M)  (compressional  attenuation  ratio  yj,ic)),  G3(1M)  (shear  attenuation  ratio  y(„is)),  GA(1M) 
(water  absorption  an),  GT(IM)  (air/water  scattering  ratio  r0,„),  and  GB(IM)  (water/sediment  scatter¬ 
ing  ratio  T!  „). 

The  values  of  the  modal  eigenfunctions  u'„(0  at  the  source  position  are  determined  as  follows. 
The  normalized  source  depth  at  range  RAN(I)  is  calculated.  Then  the  value  of  the  eigenfunction, 
XXR1,  for  this  normalized  depth  at  the  SR  range  is  determined  by  interpolating  between  stored  values 
of  the  eigenfunction  computed  for  the  SR  range  (CSN  202).  The  method  is  similar  to  that  used  for  the 
receiver  eigenfunction  which  is  given  by  Eqs.  (14)  and  (15)  and  which  is  illustrated  in  Fig.  4.  The 
interpolation  is  repeated  (CSN  203)  on  the  eigenfunction  calculated  at  the  LR  range  to  obtain  XXR2. 
Finally,  the  value  of  the  eigenfunction  at  RAN(I),  programmed  as  SM(J2,1M)  is  determined  by  a 
range-weighted  interpolation  between  XXR1  and  XXR2  at  CSN  204  where  the  interpolation  coefficient 
A  is  determined  from  Eq.  (16).  After  this  procedure  has  been  applied  to  all  the  normal  modes  which 
propagate  to  the  receiver  and  the  loss  at  RAN(I)  is  determined,  the  normalized  source  depth  at 
RAN0+ 1)  is  calculated  and  the  above  procedure  is  repeated. 

Note  that,  in  general,  XXRt  and  XXR2  will  change  as  RAN (I)  moves  between  SR  and  LR. 
Interpolation  employing  the  normalized  depth  variable  has  been  found  to  be  more  accurate  than  direct 
interpolation  using  the  depth  variable  z  to  obtain  different  normalized  depths  at  the  SR  and  LR  ranges. 

Part  8:  Transmission  Loss:  I. 

The  array  PL(J1,J2,1)  is  used  for  the  incoherent  transmission  loss  or  transmission-loss  anomaly 
(see  Part  9  for  definition  of  transmission-loss  anomaly)  and  QC(J1,J2)  and  QS(J1,J2)  are  used  for  the 
cosine  and  sine  terms,  respectively,  of  the  coherent  transmission  loss  or  transmission-loss  anomaly. 
These  arrays  are  first  set  to  zero.  Note  that  for  each  new  calculation-range  point,  they  will  initially  con¬ 
tain  all  zeros.  On  the  other  hand,  the  arrays  SE(IM) . SB(1M),  were  initialized  to  zero  outside  the 

DO  300  range  loop  (see  Part  3).  Thus  they  initially  contain  zeros  only  for  the  first  range-point  calcula¬ 
tion.  For  each  individual  mode  (they  are  subscripted  for  mode  index),  these  arrays  will  accumulate 
"pieces"  of  the  range  integrals  they  represent  (see  below  and  also  Fig.  6)  as  execution  of  the  code  con¬ 
tained  in  the  range  loop  is  repeated  for  each  new  calculation-range  point. 
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As  noted  in  the  discussion  of  Part  7,  the  average  attenuation  coefficient  A„  (see  Eq.  08))  may  be 
broken  into  a  number  of  integrals  representing  the  separate  attenuation  mechanisms.  Considering  only 
the  sediment  attenuation  ratio  y,l2),  for  example,  the  coefficient  A^  (a  component  of  A„)  is  given  by 
Eq.  (19).  Similar  equations  hold  for  the  other  attenuation  mechanisms.  The  program  approximates 
each  of  these  integrals  by  dividing  it  into  "pieces,"  each  "piece"  representing  the  integral  over  the  range 
from  RAN(l-l)  to  RAN(I);  recall  I  is  the  index  of  the  DO  300  range  loop.  The  piece 
J«t>  yj2)  (r')  dr'  is  approximated  by  the  area  of  a  rectangle  having  sides  DR  and  Gl(IM)  (see  Part 
7  and  Fig.  6).  Thus  it  is  programmed  as  G1  (1M)  *  DR.  In  the  FORTRAN  statement. 

S  KIM )  =  S  KIM )  +  G  KIM )  *  DR, 

Sl(IM)  on  the  left-hand  side  represents  J0RAN<I)  y^2*  (/)  dr'.  On  the  right-hand  side,  SKIM)  represents 
the  accumulated  value  of  the  integral  for  previous  ranges,  i.e.,  it  contains  the  approximation  for 

JRAN,.-,,  y  (2)  (r0  dr, 

Other  variables  handled  in  an  identical  manner  to  that  of  Sl(IM)  are:  SE(1M),  which  represents 
/  o  k„(r')dr'\  S2(1M),  which  represents  (r')  dr'\  S3(1M),  which  represents 

Joyi3s)(r')  dr'\  SA(IM),  which  represents  Jo  an(r')  dr' ;  ST(IM),  which  represents  J0r  r0n(/-')  dr'-, 
and  SB(IM),  which  represents  Jo  r,  „  (r')  dr' . 


Equations  (10)  and  (11)  for  the  transmission  loss  each  contain  the  term  e 
and  (5),  we  have 


r  =  «4  J0  y™  (»■')  dr'  +  es  Jg  y{nic)  ( r ')  dr'  +  e6  JQ  y,*3s)  (r')  dr' 
+  /o  S°-»  )  dr  +  Jo  ^  dr'  +  /o  a" )  dr  ' 


Using  Eqs.  (18) 


(20) 


where  we  have  inserted  e4,  «5,  and  «6  in  place  of  e2i  «jo  and  e3s,  respectively,  as  discussed  in  Part  2. 
The  term  A„  r  is  programmed  as  QQ.  We  first  encounter  the  definition 

QQ=EP4*  Sl(IM)  +  EP 5 ♦  S2(IM)  +  EP6*  S3(IM)  +  SA(1M), 


which  adds  the  first  three  terms  and  the  last  term  on  the  right-hand  side  of  Eq.  (20).  The  remaining 
two  terms  of  Eq.  (20)  are  added  to  QQ  by  the  two  FORTRAN  statements  following  the  defining  state¬ 
ment.  The  terms  S0  „  and  S,  „  may  be  expressed  in  terms  of  the  scattering  ratios  r0i„  and  T ,  respec¬ 
tively;  the  relationship  is  given  by  Eqs.  (13).  The  plane-wave  reflection  coefficients  appearing  there 
may  be  evaluated  in  terms  of  the  rats  roughnesses  of  the  boundaries.  If  we  let  a0  (S1G0  in  the  pro¬ 
gram)  be  the  rms  wave  height  and  <rx  (SIG1  in  the  program)  be  the  rms  excursion  of  the 
water/ sediment  interface,  then  Eqs.  (13)  take  the  form  (31: 


So.n 

St.„ 


2oTq 

2  W 


0> 

2 

-  k2 

c,(0) 

0,n 


O) 

2 

-  k2 

c,(//,) 

(21) 


If  the  explicit  expressions  for  r0.„  and  T i,„  given  in  Ref.  1  are  inserted  into  Eqs.  (21),  the  result  is  that 
of  Ref.  3.  The  fourth  term  on  the  right-hand  side  in  Eq.  (20)  may  now  be  written 


//j0.„  (,')*•-  2W/0' 


c,(0) 


-  kn2(r’} 


r0,„  (/■')  dr'. 
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The  FORTRAN  statement  which  includes  this  term  in  QQ  is: 

QQ  -  QQ  +  2.0«ST  (IM  )*SIG  0*SIG0*  ((6.283 1853*F/CTS  )**  2  —  WN  (IM  )*WN  (1M  )). 


The  inclusion  of  the  water/sediment  scattering  term  follows  in  a  similar  manner. 

The  term  e  X"r  is  programmed  as  Q2: 

Q2=  1.0/EXP  (QQ). 


If  A„r  >  32.25  for  any  given  mode  n,  the  term  Q2  will  not  be  included  in  the  modal  sum  of  Eqs.  (10) 
and  (11),  and  the  program  prints  out  a  flag  informing  the  user  of  this  fact.  The  reason  for  neglecting 
the  term  is  as  follows.  Suppose  that  QQ  >  32.25.  Then  Q2  <  10-14.  The  largest  value  that  Q1  can 
reasonably  be  expected  to  attain  is  of  order  104;  thus  even  for  the  largest  value  of  Ql,  Q  will  always  be 
less  than  10~10  (see  below  for  definitions  of  Ql  and  Q).  Such  a  small  term  will  not  make  a  significant 
contribution  to  the  modal  sum  appearing  in  Eqs.  (10)  and  (11),  thus  it  is  neglected.  What  this  means 
physically  is  that  at  ranges  for  which  A„r  >  32.25,  nearly  all  of  the  energy  initially  in  the  nth  mode  has 
been  removed  by  attenuation. 

Also  appearing  in  the  modal  sums  of  Eqs.  (10)  and  (11)  is  the  term  u„  (£0)  u„’  (£)/<£ ,|/2.  Except 
for  inclusion  of  a  factor  rl/J  in  the  numerator,  this  term  is  programmed  as  Ql:  Ql  “ 
RE(J1,IM)»  SM(J2,IM)»  DSQRT(RAN(I)/SE(IM)).  (See  Part  9  for  programming  of  the  removal  of 
the  factor  r'/2.) 


The  variable  Q  is  defined  as  the  product  of  Ql  and  Q2.  We  have  seen  above  that  Q  may  be 
neglected  for  a  particular  mode  if  it  is  smaller  than  10~10.  This,  in  turn,  imposes  a  restriction  on  the 
smallness  of  Q2  corresponding  to  the  largest  possible  value  of  Ql.  Alternatively,  consider  the  largest 
possible  value  of  Q2,  which  is  one.  We  may  therefore  restrict  the  absolute  value  of  Ql  to  be  greater 
than  or  equal  to  1(T10.  The  test  for  Q2  <  10~i4,  when  fulfilled,  will  remove  a  given  modal  contribu¬ 
tion  for  all  source  and  receiver  depths.  The  test  for  Ql  >  10'10  must  be  a  function  of  mode-index, 
source  depth,  and  receiver  depth.  The  corresponding  FORTRAN  statement  thus  appears  within  the 
DO  220  source  and  receiver  loops  as  well  as  within  the  DO  230  mode  loop. 

The  inclusion  of  terms  excluded  by  the  two  above-described  tests  would  not  invalidate  the 
transmission-loss  calculations.  One  reason  for  having  the  tests  is  to  save  execution  time.  However,  a 
more  important  reason  exists.  The  argument  of  the  exponential  function,  QQ,  may  become  quite  large. 
If  no  check  is  made,  a  fatal  execution  error  may  result,  due  to  machine  limitations  on  the  size  of  the 
argument  of  the  EXP  function.  Similarly,  Q2  may  become  very  small  in  absolute  value.  This  will 
occur,  for  example,  if  either  of  the  depths  Co  or  C  is  such  that  the  eigenfunction  u„  (Co)  or  (C)  is  very 
close  to  a  node  (for  a  particular  mode  order  «).  The  resulting  calculation  of  Ql  might  yield  a  number 
smaller  than  10"*,  where  x  specifies  the  dynamic  range  of  a  real  constant  (a  machine-dependent  param¬ 
eter).  This  would  cause  a  fatal  error.  Such  a  problem  has  never  been  encountered  in  years  of  using 
the  program.  On  the  other  hand,  numbers  have  been  encountered  which,  when  squared  (as  Eqs.  (10) 
and  (11)  require) ,  would  have  caused  an  execution  error  due  to  their  extreme  smallness.  The  test  on 
Ql  eliminates  the  possibility  of  such  errors. 


After  the  definition  of  Q  comes  the  FORTRAN  statement  QS(J1,J2)  =  QS(J1,J2)+Q» 
DSIN(SEOM)).  For  each  source  depth  and  each  receiver  depth  (for  which  these  arrays  are  sub¬ 
scripted),  QS(J1,J2)  will  "accumulate"  NMODE  terms  as  the  DO  230  loop  is  executed.  After  the  loop 
is  firished,  we  have: 


QS  (J  1,J  2) 


N 

I 

<7—  ) 


u„  (£0)  u'n  «)  r 


1/2 


~A_r 


1/2 


e  "  sin  <f>„. 
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< 


In  a  similar  manner,  we  have: 


0CU1.J2)  -  £ 

n*  1 


*  Un  (£0>  (£>  r1/2 


<*> 


1/2 


COS  <£„ 


and 


PL  (J  1,J  2,1 ) 


/V 


I 


K  <£o>  U'n  (£) 


,  1/2 


(22) 


The  array  PL(J1,J2,I)  will  finally  represent  the  incoherent  transmission  loss  or  transmise>on-loss  ano¬ 
maly  (see  Part  9).  As  presently  defined  by  Eq.  (22),  however,  it  represents  an  intermediate  result. 
The  multiplicative  factor  r  has  been  introduced  for  numerical  purposes  and  is  to  be  removed  in  the  final 
calculation. 


Pari  9:  Transmission  Loss:  II 

The  array  COPL(Jl,J2,I)  is  used  for  coherent  loss  calculations  and  is  defined  initially,  for  each 
source  and  receiver  depth  (via  the  DO  270  loops),  as  the  sum  of  the  squares  of  QS(J1,J2)  and 
QC(J1,J2): 


COPL  (J  1,J  2,1) 


A  U„Uo)u'(0  -A .  , 

I  - 7172 - e  s,n<*« 

n- 1  Vn 


r 


•f" 


£  u„(£o)  «„'(£)  -A„ 

I  - 7171 - * 

<1-1  9>n 


cos  4>n 


(23) 


The  phase  of  the  signal,  as  defined  in  the  "THEORY"  section,  is  coded  next. 

PHASE  (J1,J2)=ATAN2(QS  (J1,J2),QC  (J  1 , J  2))*  57.295779513. 

The  multiplicative  factor  at  the  end  converts  the  result  of  the  arctangent  function  from  radians  to 
degrees. 

The  program  will  calculate  either  the  transmission  loss  (TL)  or  the  transmission- loss  anomaly 
(TLA),  which  is  the  loss  in  addition  to  that  caused  by  cylindrical  spreading.  The  input  parameter 
ISPRD  controls  this  choice.  If  ISPRD  =  0,  then  the  TLA  is  calculated.  If  ISPRD  =  I,  then  the  TL  is 
calculated.  The  relation  between  the  two  is:  TL  =  10  log)0  r  +  TLA.  We  shall  discuss  here  only  the 
coherent  and  incoherent  TL,  which  are  given  by  Eqs.  (10)  and  (11),  respectively.  The  TLA  is  obtained 
in  a  similar  manner. 

To  get  Eq.  (10)  from  Eq.  (23),  we  divide  by  the  range  r=RAN(I),  the  water  depth  at  the  receiver 
H10,  and  the  water  depth  at  the  source  HS;  we  multiply  by  2 np}\  and  we  take  the  base-10  log  of  the 
result  and  multiply  by  (—10).  The  FORTRAN  statement  which  does  this  therefore  completes  the 
evaluation  of  Eq.  (10)  and  stores  the  final  result  in  the  array  COPL(Jl,J2,l).  To  get  Eq.  (11)  from 
Eq.  (22),  we  follow  an  identical  procedure.  The  array  PL(J1,J2,I)  contains  the  final  result.  In  writing 
Eqs.  (10)  and  (11),  the  "normal"  situation  of  a  receiver  in  the  water  layer  was  assumed.  In  this  case, 
we  have  DEPRE(J1).LE.H10  fulfilled  in  the  FORTRAN  IF  statement,  and  Eqs.  (10)  and  (11)  are  pro¬ 
grammed  with  pi  given  by  RERHOl.  As  mentioned  earlier,  however,  the  program  will  accommodate  a 
receiver  in  the  sediment.  In  this  case,  the  program  replaces  pi  by  p 2,  given  in  the  FORTRAN  by 
RERH02. 
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Part  10:  Primed  Output 

The  calculated  coherent  transmission  loss  (modulus  in  decibels,  phase  in  degrees)  and  incoherent 
transmission  loss  (in  decibels)  are  printed  out  as  functions  of  range  (DO  300),  source  depth  (DO  290), 
and  receiver  depth  (DO  280).  (When  there  is  only  one  source  and  one  receiver,  the  quantities  are 
functions  of  range  only,  and  the  output  is  on  one  line  per  range.) 

Examples  are  given  in  the  "OUTPUT"  section. 

Part  1 1:  Plotted  Output 

Plotting  of  loss  vs  range  is  performed  separately  for  each  source  depth  and  receiver  depth.  The 
plots,  which  contain  both  coherent  and  incoherent  loss,  are  executed  in  the  following  order,  which  is 
important  to  note  since  the  plots  are  unlabeled: 

DO  340  loop  for  frequency 
DO  330  loop  for  receiver  depth 
DO  330  loop  for  source  depth 
Plot  Package 

330  CALL  ORIGIN  (XLENG  +2.5,  0.) 

340  CONTINUE 


INPUT  DATA 

Explanation  of  Data  Deck  and  Notes  to  the  User 

There  are  three  groups  of  data  input  statements.  The  first  group,  containing  primarily  data 
needed  for  transmission-loss  calculations,  is  read  in  once  (in  Part  1).  (The  various  "Parts"  of  the  pro¬ 
gram  are  defined  in  the  "Step-By-Step  Analysis"  section.)  The  second  group,  containing  primarily 
environmental  data  needed  for  modal  calculations  at  the  receiver  site,  is  read  in  once  for  each  fre¬ 
quency  (in  Part  3).  The  third  group  is  similar  to  the  second  and  is  read  in  (only  for  a  range-dependent 
calculation)  NRBUF-1  times  for  each  frequency  (in  Part  6).  Each  of  the  second  and  third  groups  may 
differ  from  the  others  by  any  arbitrary  combination  of  the  environmental  parameters,  subject  to  the 
constraint  that  the  fluid-basement  and  solid-basement  models  may  not  both  be  used  for  a  given  track. 
If  NMFREQ  is  greater  than  one,  the  procedure  of  reading  the  second  and  third  groups  is  repeated. 

The  spacing  of  the  user-supplied  environmental  data  profiles  for  the  range-dependent  model  is 
usually  dictated  by  the  bathymetry  of  the  acoustic  propagation  path,  and  to  a  lesser  degree,  by  changes 
in  the  sound  speed  profile  with  range.  Since  water  depth  at  a  calculation  range  site  is  approximated  by  a 
linear  interpolation  between  the  depths  at  two  user-supplied  profiles,  the  user  should  approximate  the 
known  bathymetry  by  linear  segments,  supplying  input  data  at  the  ranges  where  these  segments  meet. 
Additional  profiles  may  be  supplied  at  other  ranges  in  order  to  take  account  of  the  range-dependence  of 
other  environmental  parameters. 

There  is  another  reason  for  inserting  extra  profiles  along  the  track.  If  the  number  of  mode;  which 
the  environment  will  support  decreases  rapidly  with  increasing  range,  as  will  be  the  case  for  a  water 
depth,  H j,  which  decreases  rapidly,  then  the  calculated  loss  may  change  discontinuously  at  an  input 
range  where  there  are  many  fewer  modes  than  at  the  previous  input  range.  This  is  due  to  conservation 
of  mode-index,  implemented  in  the  program  by  excluding  those  higher  order  modes  which  are  not 
present  at  both  of  the  input  environments  between  which  range-interpolations  are  to  be  carried  out  on 
the  eigenfunctions.  Thus,  for  example,  if  three  input  profiles  A,  B,  and  C  support  80.  50,  and  20 
modes,  respectively,  calculations  at  range  points  between  A  and  B  will  use  50  modes,  and  those 
between  B  and  C  will  use  20  modes.  At  ranges  greater  than  B,  some  of  the  (excluded)  higher  order 
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modes  may  initially  actually  yield  an  important  contribution  which  decreases  rapidly  thereafter.  By 
exclusion  of  all  these  modes  starting  at  range  B,  a  small  discontinuity  in  the  loss  may  be  caused  at  range 
B.  The  problem  may  be  avoided,  however,  by  supplying  additional  profiles,  perhaps  one  between  A 
and  B  which  supports  approximately  65  modes  and  one  between  B  and  C  which  supports  approximately 
35  modes.  (If  the  discontinuity  is  reduced,  but  not  eliminated,  more  supplemental  profiles  may  be 
necessary.) 

The  abovi-described  procedure  will  work  if  enough  intermediate  profiles  are  used;  however,  this 
may  substantially  increase  the  execution  time  of  the  program.  We  shall  briefly  discuss  another  method 
of  circumventing  'he  problem  which  is  particularly  suited  to  those  tracks  along  which  the  deep  sediment 
structure  is  not  important.  Reference  to  Eq.  (12)  shows  that  the  total  number  of  modes  which  exists  at 
a  given  range  is  dt.'endent,  among  other  things,  on  H  =  H,  +  H2,  and  not  on  the  water  depth  Hi 
alone.  (It  was  tacitly  >ssumed  in  the  previous  paragraph  that  H2  is  nearly  constant  with  range  so  the 
change  in  the  number  of  allowed  modes  was  due  to  the  decrease  of  H\.)  If  deep  sediment  structure  is 
unimportant,  the  sediment  thickness  H2  may  be  assigned  small  values  at  the  smaller  ranges  and  large 
values  at  the  larger  ranges,  thereby  keeping  H  and  consequently  the  total  number  of  modes  nearly  con¬ 
stant  with  changing  range.  This  procedure  allows  the  higher  order  modes  to  accumulate  large  attenua¬ 
tion  coefficients  before  their  contributions  are  dropped  from  the  modal  sum  in  Eq.  (22). 

In  experiments  the  source  depth  sometimes  changes  at  one  or  more  points  along  the  tow  track. 
The  user  of  PROGRAM  MOATL  can  incorporate  source  depth  changes  along  the  track  by  inserting  the 
appropriate  IF  statement(s)  between  CSNs  193  and  194: 

IF(RAN(I).GT.  . )  DEPSO(l)  =  ...  (etc). 

In  situations  in  which  the  water  depth  over  part  of  the  acoustic  propagation  path  is  less  than  the  source 
depth,  the  change  of  source  depth  is  necessary  since  the  program  requires  the  source  to  be  in  the  water 
layer.  This  change  is  necessary  even  if  the  source  was  not  towed  over  that  part  of  the  path  in  the 
experiment  whose  results  are  to  be  modeled,  since  the  program  assumes  that  the  source  is  towed  from 
the  range  origin  in  order  to  evaluate  Eqs.  (18),  (19),  and  (20)  accurately.  We  note  for  completeness 
that  for  the  range-dependent  or  range-independent  model,  the  source  depth  may  be  changed  at  any 
range  for  any  other  reason  and  all  the  loss  calculations  will  be  correct. 

There  are  two  other  notable  variations  that  the  program  will  treat.  First,  although  the  program  is 
set  up  to  model  a  three-layer  half-space,  a  two-layer  model  may  be  practically  realized  when  the  base¬ 
ment  rock  interfaces  with  the  water  (SOLID)  or  the  sediment  is  so  thick  that  deep  sediment  structure  is 
not  important  (FLUID).  Such  an  example  for  the  fluid  basement  is  given  in  the  "OUTPUT'  section. 
The  second  variation  is  that,  although  the  code  as  programmed  assumes  the  calculation  ranges  to  be 
equally  spaced,  minor  code  alterations  will  allow  unequally  spaced  ranges.  Details  concerning  each  of 
these  two  variations,  along  with  additional  information,  may  be  found  in  the  section  "General 
Remarks." 

Listed  below  are  the  required  input  data.  Each  line  (and  an  indented  continuation)  corresponds  to 
a  single  data  card  (or  a  single  record  of  80  characters  if  input  is  not  by  use  of  cards).  If  there  is  more 
than  one  of  each  type  of  card,  the  variable  specifying  the  number  of  cards  is  printed  to  the  left  and  is 
underlined  All  of  such  cards  appear  together  in  the  data  deck.  (By  "NDSO  -t-8”  it  is  meant,  for  exam¬ 
ple,  that  if  there  are  to  be  eleven  source  depths,  there  will  be  two  data  cards,  eight  values  on  the  first 
and  three  on  the  second.)  The  three  input  groups  described  in  the  first  paragraph  of  this  section  are 
separated  by  brackets  If  there  is  more  than  one  set  of  a  group  of  cards,  the  variable  specifying  the 
number  of  sets  is  printed  to  the  left  and  is  underlined.  These  sets  are  clustered  in  the  data  deck.  The 
FORMAT  to  be  used  for  each  card  is  specified  in  parentheses  at  the  end  of  each  line. 
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NMFREQ 


TITLE(I)  (20A4) 

1SPRD,NMFREQ,NRBUF,NRCALC,RMAX,EP1,EP2, 
EP3,SIG0,SIG1  (4I5,F10.3,SF10.7) 

NDSO,NDRE  (215) 

HI  (15) 

lPLOT,DBMIN,DBMAX,DY,DX  (I5,2F10.3,2  F10.7) 
NDSO±8  -  DEPSO(l),  1-1,8  (8F10.3) 

NDRE-r8  -  DEPRE(l),  1-1,8  (8F10.3) 

TITLE(I)  (20A4) 

F  (F10.3) 

MDPRNT,1NC1,1NC2,RH01,RH02,RH03,H1 1  ,H2 
(315,5F10.3) 

EPSLN,COMP,SHEAR,RANGEl  ,LI1  ,L12,ND1  ,ND2 
(F10.8.3F1 0.3,415) 

NDi-Zl(l),Cl(l)  (2F10.3) 

ND2  — Z2(I),C2(I)  (2F10.3) 


NRBUF-1 


MDPRNT,INC1,INC2,RH01,RH02,RH03,H12,H2 

(3!5,5F10.3) 

EPSLN, COMP, SHE  AR,RANGE2,L11,LI2,ND1,ND2 
(F10.8,3F10.3,4I5) 

NDi -Z1(1),C1(I)  (2F10.3) 

ND2  —  Z2(I),C2(I)  (2F10.3) 
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Description  of  Input  Data 

We  now  describe  the  function  and  use  of  each  of  the  input  variables.  They  are  arranged  here  in 
the  same  order  as  they  are  read  in  by  the  program.  For  determining  input  values  further  information  is 
given  in  the  previous  section  and  the  "General  Remarks"  section. 

No  specific  units  (m  or  km,  for  example)  are  required;  however,  it  is  necessary  to  see  that  the 
units  chosen  are  consistent  for  all  the  input  variables.  (The  only  two  exceptions  are  the  density  of  the 
water,  which  must  always  be  unity  (see  below),  and  DX  (km/in.)  and  RMAX  (m)  when  the  plot  pack¬ 
age  is  desired.)  We  assume  below  the  m-k-s  system  for  convenience,  which  serves  also  as  a  useful 
example. 

TITLE(I)  —  An  array  containing  any  alphanumeric  label  (80  characters  maximum). 

ISPRD  —  If  0,  the  loss  anomaly  (no  cylindrical  spreading)  will  be  calculated.  If  1,  the  loss  (cylindrical 
spreading  included)  will  be  calculated.  (See  discussion  in  the  section  "Step-By-Step  Analysis- 
Part  9".) 

NMFREQ  —  The  number  of  source  frequencies  (or  alternatively,  the  number  of  different  environmen¬ 
tal  cases  for  a  given  frequency)  for  which  the  calculations  are  desired.  Should  be  set  to  unity 
when  unequally  spaced  calculation  ranges  are  used. 

NRBUF  —  The  number  of  environmental  data  profiles  to  be  input.  If  1,  the  environment  is  range- 
independent.  If  greater  than  1,  the  environment  is  range-dependent. 

NRCALC  —  The  number  of  range  points  at  which  loss  calculations  are  desired  (limited  by  PARAME¬ 
TER  statement).  In  general,  these  will  be  equally  spaced  out  to  the  maximum  range  specified. 
(See  the  "General  Remarks"  section  for  further  discussion  and  limitations.) 

RMAX  —  Maximum  range  (m)  at  which  loss  calculations  are  desired. 

EP1  —  Sediment  layer  plane-wave  absorption  coefficient  (dB/Hz-m). 

EP2  —  Basement  compressional  plane-wave  absorption  coefficient  (dB/Hz-m). 

EP3  —  Basement  shear  plane-wave  absorption  coefficient  (dB/Hz-m).  Set  equal  to  0  for  a  FLUID 
model. 

SIGO  —  The  rms  wave  height  (m)  used  for  calculating  attenuation  due  to  air/water  interface  scattering. 

SIGI  —  The  rms  bottom  roughness  (m)  used  for  calculating  attenuation  due  to  water/sediment  inter¬ 
face  scattering. 

NDSO  —  Number  of  source  depths.  (Limited  by  PARAMETER  statement.) 

NDRE  —  Number  of  receiver  depths.  (Limited  by  PARAMETER  statement.) 

Ill  —  If  0,  a  FLUID  basement  is  assumed.  If  1,  a  SOLID  basement  is  assumed.  (The  type  of  base¬ 
ment  is  fixed  by  III  for  the  entire  track.) 

IPLOT  —  If  0,  no  plots  are  included  in  the  output.  If  1,  plotting  is  executed. 

DBMIN  —  Lower  bound  for  transmission  loss  (dB)  on  the  plotted  graph’s  loss  axis. 
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DBMAX  —  Upper  bound  for  transmission  loss  (dB)  on  the  plotted  graph's  loss  axis. 

DY  —  Number  of  decibels  per  inch  determining  the  scale  of  the  plotted  graph’s  loss  axis. 

DX  —  Determines  the  scale  of  the  plotted  graph’s  range  axis.  If  the  units  used  for  range  are  meters, 
DX  should  specify  the  number  of  kilometers  per  inch  on  the  axis. 

DEPSO(l)  —  An  array  storing  the  source  depths  (m).  (Each  of. these  values  must  be  <  Hf.) 

DEPRE(l)  —  An  array  storing  the  receiver  depths  (m).  (Each  of  these  values  must  be  <  H:  +  H 2.) 

TlTLE(l)  —  Another  arbitrary  alphanumeric  label  (80  characters  maximum). 

F  —  Source  frequency  (Hz). 

MDPRNT  —  If  0,  none  of  the  calculated  mode-amplitude  functions  will  be  printed  out.  Otherwise, 
they  will  be  printed  out  in  accordance  with  the  following  two  inputs. 

INC1  —  A  variable  containing  a  value  one  greater  than  the  number  of  depth-increment  values  to  be 
skipped  between  printed-out  values  of  the  mode-amplitude  functions,  in  the  water  layer.  Note: 
first  determine  LI  1 ;  then  set  1NC1  to  give  the  desired  number  of  printed  out  values.  As  an 
example,  if  Lll  =  300  and  INC1  =  6,  then  50  values  will  be  printed  out. 

INC2  —  Same  as  INC1,  but  for  the  sediment  layer.  The  number  of  output  values  is  based  on  LI2. 

RHOl  —  The  variable  RHOl^l.O  always.  It  is  the  density  to  be  used  for  water  regardless  of  the  units 
used  for  the  rest  of  the  data. 

RH02  —  Ratio  of  the  density  of  the  sediment  layer  to  that  of  water. 

RH03  —  Ratio  cf  the  density  of  the  basement  to  that  of  water. 

Hll  —  The  thickness  (m)  of  the  water  layer  at  the  receiver. 

H12  —  The  thickness  (m)  of  the  water  layer  at  any  nonzero  range  where  the  environmental  data  are 
read  in. 

H2  —  The  thickness  (m)  of  the  sediment  layer.  (Same  variable  name  used  at  all  ranges.) 

EPSLN  —  The  variable  EPSLN  =  0.0001  always.  It  is  the  criterion  for  the  accuracy  of  the  calculated 
eigenvalues  (modal  wave  numbers),  i.e.,  the  amount  by  which  the  air/water  surface  value  of 
the  normalized  eigenfunction  may  differ  from  the  pressure-release  boundary-condition  require¬ 
ment  of  being  identically  zero. 

COMP  —  Compressional  velocity  (m/s)  in  the  basement. 

SHEAR  —  Shear  velocity  (m/s)  in  the  basement  for  the  SOLID  model;  it  must  exceed  the  minimum 
sound  speed  found  in  the  water  and  sediment  layers.  However,  for  the  FLUID  mode),  set 
SHEAR  =  0. 

RANGE1  —  The  variable  RANGE1  =  0  always.  It  is  the  receiver  range,  which  is  zero  by  defintion. 

RANGE2  —  The  range  (m)  from  the  fixed  receiver  of  the  present  environmental  profile  being  read  in. 
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Lll  —  The  number  of  incremental  steps  into  which  the  water  layer  is  to  be  divided  for  the  calculations. 
(See  "Notes"  below.) 

LI2  —  The  number  of  incremental  steps  into  which  the  sediment  layer  is  to  be  divided  for  the  calcula¬ 
tions.  (See  "Notes"  below.) 

Notes:  The  numbers  of  layers  Lll  and  L12  help  to  determine  the  accuracy  of  the  eigenvalues  and 
mode-amplitude  functions.  The  quantity  Lll  should  be  set  equal  to  about  ten  times  the  number  of 
modes  expected;  L12  should  be  set  to  yield  approximately  the  same  spacing.  To  estimate  N,  the 
total  number  of  modes,  see  Eq.  (12).  Dimensioning  of  arrays  requires  that  N  be  less  than  or  equal 
to  150.  In  no  event  should  either  Lll  or  LI2  be  less  than  four,  nor  should  their  sum  exceed  1200. 
They  must  each  be  even-valued  to  be  consistent  with  the  Simpson’s  rule  integration  method  used 
in  the  subroutines. 

ND1  —  The  number  of  sound-speed  profile  depths  to  be  read  in  for  the  water  layer  (not  to  exceed 
150). 

ND2  —  The  number  of  sound-speed  profile  depths  to  be  read  in  for  the  sediment  layer  (not  to  exceed 
150). 

Z1  (1) ,  CKI)  —  The  arrays  for  the  depth  (m)  and  sound  speed  (m/s),  respectively,  profile  values  of  the 
water  layer. 

Z2(I),  C2(I)  —  The  arrays  for  the  depth  (m)  and  sound  speed  (m/s),  respectively,  profile  values  of  the 
sediment  layer. 

Notes.  The  following  conditions  must  be  met:  Zl(l)  =  0;  Zl(NDl)  =  Z2(l)  =  Hll  (or  H12,  if 
range  r  ^  0);  Z2(ND2)  —  Hll  +  H2  (or  H12  +  H2).  Since  the  program  interpolates  linearly  for 
sound  speeds  between  those  given,  it  is  sufficient  to  supply  only  the  two  bounding  depths  when 
there  are  three  or  more  consecutive  depth  points  having  sound  speed  a  linear  function  of  depth. 

For  the  user  who  knows  only  the  type  of  material  comprising  the  sediment,  we  include  Table  2,  as 
a  guide  in  determining  input  values.  The  density  of  a  given  sediment  type  is  suitable  for  the  input  vari¬ 
able  RH02  as  given.  The  velocity  ratio,  when  multiplied  by  the  value  of  Cl(NDl),  yields  the  value  to 
be  used  for  C2(l).  If  the  attenuation  coefficient  (dB/m-kHz)  is  to  be  used  for  EP1  (supplied  in  dB/m- 
Hz),  it  must  first  be  multiplied  by  10~3. 


Table  2°  —  Sediment  Layer  Parameters 


Sediment 

Type 

Density 

(g/cm3) 

Porosity 

(%) 

Velocity 

Ratio 

Atten.  Coeff. 
(dB/m-kHz) 

Coarse  sand 

2.034 

38.6 

1.201 

0.47 

Fine  sand 

1.957 

44.8 

1.147 

0.51 

Very  fine  sand 

1.866 

49.8 

1.111 

0.68 

Silty  sand 

1.806 

53.8 

1.091 

0.69 

Sandy  silt 

1.787 

52.5 

1.088 

0.76 

Silt 

1.767 

54.2 

1.062 

0.68 

Sand-silt-clay 

1.583 

67.2 

1.033 

0.11 

Clayey  silt 

1.469 

72.6 

1.011 

0.08 

Silty  clay 

1.421 

75.9 

0.994 

0.07 

“Compiled  by  Anthony  1.  Eller  and  I  rank  Ingenito  Trom  data  given  in  Refs  9- 
12  —  pi \\aw  communication 
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OUTPUT 

Most  of  the  input  variables  are  also  printed  out,  via  WRITE  statements  that  follow  the 
corresponding  READ  statements  in  the  code.  This  helps  in  error-checking  and  also  serves  to  label  the 
printout  uniquely.  The  first  page  of  printed  output  contains  all  of  the  input  variables  of  group  one  (the 
three  groups  of  input  variables  are  discussed  in  the  "INPUT  DATA"  section)  except  for  the  variable  III. 
The  value  of  III  is  reflected,  however,  in  the  inclusion  or  omission  (on  the  second  printed  page)  of  the 
shear  and  Rayleigh  velocities.  Apart  from  the  latter,  which  is  a  calculated  quantity,  the  second  page  of 
output  contains  only  the  values  of  variables  of  the  second  input  group. 

Beginning  on  the  third  page  are  the  calculated  modal  properties  of  each  of  the  normal  modes  sup¬ 
ported  by  the  environment  given  on  the  second  page.  The  FORTRAN  WRITE  commands  for  this  out¬ 
put  are  located  near  the  end  of  SUBROUTINE  SOLID  (or  FLUID).  For  each  mode,  the  mode-order 
and  phase  velocity  are  first  printed  out.  Also  printed  is  the  number  of  iterate  solutions,  i.e.,  the 
number  of  times  that  SUBROUTINE  HALFS  (or  HALFF)  called  SUBROUTINE  1TRTS  (or  1TRTF)  in 
order  to  converge  to  the  correct  eigenvalue.  Also  printed  on  the  same  line  is  the  number  of  times  that 
the  eigenfunction  had  to  be  scaled  down,  according  to  the  technique  described  in  Appendix  A.  The  last 
two  lines  for  a  given  mode  contain,  from  left  to  right,  the  label  and  value  for  each  of  the  following 
quantities:  the  wave  number  (eigenvalue  kn),  the  water  absorption  a„  the  sediment-layer  attenuation 
ratio  yi2),  the  basement  compressional  attenuation  ratio  y,J3c),  the  basement  shear  attenuation  ratio 
y ’  (included  only  for  the  solid-basement  model),  the  air/water  scattering  attenuation  ratio  r0  n,  and 
the  water/sediment  scattering  attenuation  ratio  F|  „. 

In  addition  to  the  modal  output  described  in  the  last  paragraph,  there  are  several  statements  which 
are  conditionally  printed  for  each  mode.  The  statement  "UPPER  AMPLITUDES  MATCHED  FOR 
THIS  MODE  STARTING  AT  NORMALIZED  DEPTH  =  ”  will  be  printed  out  (along  with  the 
appropriate  value  of  the  normalized  depth  zm)  if  the  eigenfunction  u„(z)  has  been  calculated  for  z<zm 
according  to  the  procedure  described  in  Appendix  A.  The  remaining  conditional  statements  are 
"LAYER  2  ATTEN  RATIO  =  DEFAULT  ZERO"  and  similar  statements  for  the  compressional  and/or 
shear  attenuation  ratios.  These  statements  are  printed  out  only  if  during  their  calculation  they  were 
determined  to  be  small  enough  to  be  set  to  zero. 

After  the  modal  parameters  and  flags  have  been  printed  out  for  each  mode,  and  if  MDPRNT  = 
1,  the  eigenfunctions  are  printed  out.  (If  MDPRNT  «=  0,  they  are  not  printed  out.)  They  are  printed 
out  in  columnar  form,  twelve  to  a  page,  with  depth  increasing  down  the  page.  (The  first  column  con¬ 
tains  the  values  of  the  normalized  depth  for  each  line.) 

If  the  problem  is  a  range-dependent  one,  the  second  environment  {i.e.,  the  one  at  the  first 
nonzero  input  range)  and  the  normal-mode  parameters  for  it  are  next  printed  out  in  a  format  identical 
to  the  first  (receiver)  environment. 

Next  to  be  printed  out  are  the  transmission  loss  calculations  for  the  output  ranges  which  fall 
between  the  first  and  second  environments.  (For  a  range-independent  problem,  there  is  no  second 
environment,  and  all  the  transmission-loss  calculations  follow  the  output  for  the  receiver  environment.) 
If  there  is  a  third  input  range,  its  environmental  and  modal  parameters  are  printed  next,  followed  by 
loss  calculations  at  output  ranges  between  the  second  and  third  input  ranges.  This  procedure  is 
repeated  until  the  output  ranges  have  been  exhausted. 

For  those  pages  containing  the  calculated  loss,  the  coherent  transmission  loss  is  printed  on  the  left 
sides  of  the  pages;  the  incoherent  transmission  loss  is  printed  on  the  right  sides  of  the  pages.  The  for¬ 
mal  for  each  is  as  follows.  The  calculations  for  each  output  range  are  grouped  together.  Following  the 
line  on  which  the  value  of  this  range  is  printed  is  a  "group”  of  output  for  each  source  depth,  which  is 
also  printed  out.  This  "group"  will  be  only  one  line  if  there  are  five  or  fewer  receiver  depths.  (If  there 
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are  between  six  and  ten  receiver  depths,  there  will  be  two  lines,  etc.)  Each  line  contains  the 
transmission  loss  for  up  to  five  values  of  the  receiver  depth.  If  there  are  fewer  than  an  integral  multi¬ 
ple  of  five  receiver  depths,  the  line  is  "filled  out"  with  indeterminate  form  data  O.e.,  a  string  of  the 
letter  "D.  The  procedure  for  filling  variables  with  indeterminate  form  data  prior  to  execution  is 
machine-dependent  and  is  not  described  here.  The  transmission  loss  is  not  labeled  with  receiver  depth, 
but  the  depths  are  in  order  of  decreasing  depth,  as  printed  out  on  the  first  page  of  output. 

We  now  include,  on  the  following  pages,  the  computer  output  (described  in  general  above)  for 
two  dissimilar  test  cases.  The  differences,  summarized  in  Table  3,  illustrate  many  of  the  options  avail¬ 
able  to  the  user. 


Table  3  —  Summary  of  Differences  Between  Test 
Cases  1  and  2 


Item 

Test  Case  1  (for  Fig.  2) 

Test  Case  2 

Type  of  model 

Range-independent  environment 

Range-dependent  environment 

Type  of  environment 

2-layer  fluid-basement 

3-layer  solid-basement 

Surface  and 

No 

Yes 

bottom  scattering 

No.  of  output 

200 

10 

ranges 

Plots 

Yes 

No 

No.  of  source  depths 

1 

2 

No.  of  receiver  depths 

3 

2 

Receiver  in  the 

No 

Yes 

sediment  layer 

Sound  speed  profile 

in  water  and 

Isovelocity 

Depth-dependent 

sediment  layers 

Frequency 

20  Hz 

100  Hz 

No.  of  modes 

6 

19,  22,  24  (3  environments) 

Eigenfunctions  printed  out 

Yes 

No 

Note  that  test  case  1  calls  for  plots  of  transmission  loss.  These  plots  have  already  been  presented 
as  Fig.  2  (see  the  "THEORY"  section).  Note  also  that  test  case  1  asks  for  transmission-loss  calculations 
at  200  range  points.  We  include  only  the  first  and  last  page  of  calculated  results  here.  The  23  interven¬ 
ing  pages  of  output  have  been  deleted. 
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Appendix  A 

NORMAL-MODE  SUBROUTINES 


Detailed  documentation  [Al]  is  available  for  the  normal-mode  calculations.  Although  the  com¬ 
plete  model  theory  and  its  FORTRAN  coding  may  be  found  in  Ref.  Al,  we  present  here  an  outline  of 
the  salient  features,  with  particular  emphasis  on  the  major  changes  and  improvements  made  to  the  rou¬ 
tines  since  the  publication  of  Ref.  Al. 


For  a  shear-supporting  solid  basement,  PROGRAM  MOATL  calls  subroutines  SOLID,  HALFS, 
and  ITRTS  for  the  normal-mode  calculations.  (The  fluid-basement  case  is  similar,  but  less  complicated, 
and  will  not  be  discussed  here.)  The  three  subroutines  used  in  the  solid-basement  case  correspond, 
respectively,  to  PROGRAM  SOLID,  SUBROUTINE  HALF,  and  SUBROUTINE  ITERATE  of  Ref.  Al. 


The  normal  modes  of  a  given  environment  are  the  eigenfunctions  of  Eq.  (Al).  Application  of 
appropriate  boundary  conditions  will  yield  a  set  of  N  discrete  modes,  with  corresponding  modal  wave 

numbers,  or  eigenvalues  k„  (n  =  1 . N).  Due  to  the  appearance  of  the  depth-dependent  r(£), 

Eq.  (Al)  must  be  replaced  by  a  finite-difference  equivalent  for  numerical  solution;  the  water  and  sedi¬ 
ment  layers  must  therefore  be  divided  into  a  number  of  intervals,  or  incremental  steps,  and  these  are 
user-specified  through  the  input  variables  LI1  and  LI2. 


d\ 

di! 


Hi 


c(0 


-  k 2 


u„  -  0. 


(Al) 


Normal-mode  calculations  begin  in  SOLID  with  a  determination  of  the  maximum  number  of 
modes.  The  calculation  of  modal  parameters  follows.  The  mode  order  is  designated  and  upper  and 
lower  bounds  are  set  for  the  eigenvalue.  HALFS  is  then  called  to  repeatedly  shorten  the  length  of  this 
interval,  thus  converging  on  the  correct  eigenvalue.  For  each  new  trial  eigenvalue  obtained,  a  call  to 
ITRTS  is  made,  where  a  set  of  "three-point  extrapolation"  equations  (the  finite-difference  equivalent  of 
Eq.  (Al))  is  used  to  generate,  from  bottom  to  surface,  the  corresponding  trial  eigenfunction.  Depend¬ 
ing  on  the  value  of  the  trial  eigenfunction  at  the  air/water  surface,  the  trial  eigenvalue  becomes  either 
the  new  left-  or  new  right-hand  bound  of  the  now-smaller  eigenvalue  interval.  Generally  ten  or  more 
calls  are  made  to  ITRTS  before  an  acceptable  eigenvalue  is  found. 

The  correct  eigenvalue  has  an  eigenfunction  which  fulfills  the  boundary  condition  of  being  zero  at 
the  air/water  (pressure-release)  surface.  Alternatively,  HALFS  requires  that  either  of  two  criteria  are 
met:  (1)  the  surface  value  of  the  trial  eigenfunction  is  less  than  EPLON  and  the  difference  between  the 
present  trial  eigenvalue  and  either  the  left  or  the  right  bound  of  the  present  interval  is  less  than  10~12; 
or  (2)  the  difference  between  the  trial  eigenvalue  and  either  bound  of  the  interval  is  less  than  PPRCN. 
The  variable  PPRCN  is  defined  as  10-PRCSN,  where  PRCSN  is  one  less  than  the  approximate  number  of 
decimal  digits  associated  with  the  mantissa  of  a  DOUBLE  PRECISION  floating-point  number.  What 
this  means  is  that  for  criterion  (2),  the  eigenvalue  is  correct  to  within  the  limits  of  the  machine.  The 
machine-dependent  PRCSN  is  adjusted  by  a  PARAMETER  statement  in  MOATL. 

For  criterion  (1),  EPLON  is  defined  as  the  product  of  EPSLN  (a  preselected  small  number)  and 
UMAX  (the  maximum  value  that  the  unnormalized  trial  eigenfunction  takes  on).  When  the  correct 
eigenfunction  is  determined,  control  returns  to  SOLID,  where  the  eigenfunction  is  normalized  (using 
the  Simpson’s  rule  numerical  integral  technique)  according  to  the  solid-basement  equivalent  of  Eq.(A2). 
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The  value  of  EPSLN  is  usually  selected  on  the  assumption  that  the  maximum  value  of  the  eigenfunc¬ 
tion  is  of  the  order  of  magnitude  one.  Since  the  surface  value  is  checked  in  HALFS  using  the  unnor¬ 
malized  trial  eigenfunction,  EPSLN  is  redefined  as  EPLON,  as  stated  above.  In  this  way,  the  normal¬ 
ized  eigenfunction  will  have  a  surface  value  less  than  the  input  variable  EPSLN,  if  criterion  (1)  is  met. 

/o“p(£>  «„2({)  dt,  =  1.  (A2) 

In  this  connection,  we  mention  another  major  change  in  the  programming.  It  was  found  in  some 
cases  that  the  unnormalized  eigenfunction  generated  in  1TRTS  took  on  very  large  values,  sometimes 
larger  than  PMGTD,  which  is  defined  as  ioMGNTD.  (The  variable  MGNTD  depends  on  one-half  the 
dynamic  range  of  a  real  constant  and,  being  machine-dependent,  is  defined  for  convenience  by  a 
PARAMETER  statement  in  MOATL.)  In  the  course  of  normalization  of  the  eigenfunction,  certain 
single-precision  variables  are  set  equal  to  the  sum  of  the  squares  of  unnormalized  eigenfunction  values. 
Thus  a  value  larger  than  PMGTD,  when  squared,  will  cause  a  floating  point  overflow  in  the  machine, 
with  subsequent  termination  of  execution  and/or  the  generation  of  erroneous  results.  This  problem 
was  eliminated  in  the  following  way.  During  generation  of  the  trial  eigenfunction  in  1TRTS,  whenever 
a  value  larger  than  PMGTD  is  encountered,  the  entire  eigenfunction  thus  far  generated  is  scaled  down 
by  dividing  the  value  at  each  incremental  depth  by  PMGTD.  If  a  given  value  is  so  small  that  such  a 
division  would  cause  a  floating  point  underflow,  that  value  is  redefined  to  be  zero.  Then  the  generation 
of  the  rest  of  the  eigenfunction  resumes.  Whenever  such  a  scaling  of  the  eigenfunction  takes  place,  a 
counter,  1FG,  is  increased  by  one.  The  value  of  this  counter  is  passed  back  to  HALFS  along  with  the 
trial  eigenvalue  and  eigenfunction.  In  HALFS,  the  eigenvalue  interval  redefinition  (described  above) 
depends  on  the  surface  values  of  the  trial  eigenfunctions,  so  the  scaling  information  must  be  taken  into 
account. 

There  is  another  important  change  in  the  programming  of  the  eigenfunction  generation  to  be  dis¬ 
cussed.  It  is  stated  in  Ref.  A1  that  for  certain  downward-refracting  sound  speed  profiles,  the  lowest 
order  eigenfunctions  generated  may  exhibit  large  amplitude  fluctuations  near  the  surface.  Previously, 
the  eigenfunction  was  zeroed  (the  values  at  each  of  the  incremental  depths  were  set  to  zero)  from  this 
point  up  to  the  air/water  surface.  Recently,  a  method  of  calculating  the  actual,  albeit  small,  values  of 
the  function  has  been  implemented  to  replace  the  zeroing  procedure. 

The  phase  velocity,  cp„  =  a>/k„,  increases  monotonically  with  mode-order  n.  Consider  a  sound- 
speed  profile  which  is  strongly  downward  refracting  near  the  surface,  such  as  the  deep  water  profile 
shown  in  Fig.  Al.  It  may  happen  that  the  phase  velocities  of  the  lowest  order  eigenfunctions,  although 
larger  than  cmm  as  they  must  be,  are  smaller  than  the  surface  value  of  the  sound  speed  profile.  (In  Fig. 
Al,  this  is  true  of  the  nth  mode,  but  not  the  (n  +  l)th  mode.)  Thus  for  some  depth  value  zm,  we  have 

kn  >  —7-7  for  0  ^  z  <  zm  . 

The  result  is  that  solutions  of  Eq.  (A  I)  for  z  <  zm  are  exponential  rather  than  sinusoidal,  and  further 
iteration  of  the  solution  may  become  unstable. 

The  program  will  now  calculate  the  values  of  the  eigenfunctions  for  z  <  zm  (rather  than  setting 
them  to  zero).  The  technique,  for  a  given  mode,  is  simply  to  start  the  iteration  again,  this  time  at  the 
air/water  surface,  and  generate  the  eigenfunction  down  to  zm.  The  two  functions  are  then  matched  at 
z„  to  obtain  the  entire  scaled  eigenfunction. 

To  discuss  this  technique,  it  is  convenient  to  use  the  symbols  and  notation  of  Ref.  Al,  to  which 
the  reader  is  referred  for  appropriate  definitions. 

The  flag  IZERO,  which  in  the  past  controlled  the  zeroing  procedure,  now  controls  the  matching 
procedure.  It  is  set  to  one  for  the  final  call  to  SUBROUTINE  1TRTS,  which  is  the  call  that  generates 
the  acceptable  eigenfunction.  Note  that  when  this  final  call  is  made  the  eigenvalue  k„  has  already  been 
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the  corresponding  modal  phase  velocity.  For  modes  of  order  >  n,  zm  —  0. 


determined.  The  matching  procedure  is  not  used  (or  needed)  in  the  eigenvalue  search  procedure. 
When  IZERO  =  1,  the  DO  6  loop  is  executed,  which  searches  the  water  layer  for  zm.  It  may  also  hap¬ 
pen  that  kn  >  (ro/ c i  (z)J  everywhere  in  the  water  layer.  In  this  case,  the  DO  8  loop  searches  the  sedi¬ 
ment  layer  for  zm  ( zm  is  always  less  than  H\  +  H2).  The  variable  MATCH  stores  the  number  of  the 
incremental  layer  corresponding  to  zm.  Then  ITRTS  proceeds,  as  usual,  to  generate  the  eigenfunction 
from  bottom  up  to  zm.  Then  control  transfers  to  statement  number  131,  where  we  begin  a  similar  cal¬ 
culation  from  the  surface  downward.  The  flag  IZERO  is  set  to  two;  this  instructs  SOLID  to  print  a 
message  informing  the  user  that  the  present  mode  has  been  "matched"  at  depth  zm.  Next,  to 
accomodate  the  pressure-release  boundary  condition,  we  set  Zno>(0)  =  0.  We  estimate  the  value  at  the 
next  incremental  depth  down  from  the  surface  by  using  the  Taylor  polynomial  of  second  order, 
expanded  about  zero: 


Z„(1)(/,) 


Z„(l’(0)  +  f, 


dZ, 


(i) 


dz 


+ 


t±  d*zr 

2!  dz2 


Until  normalization,  the  slope  of  the  eigenfunction  is  arbitrary.  To  avoid  the  problem  of  "exponential 

dZ„(1> 

runaway,"  we  therefore  choose  the  very  small  value  of  10~'°  for 


27<1> 


dlZ, 


dz1 


dz 


We  can  solve  Eq.  (Al)  for 


,  and  the  result  is  zero.  Thus  we  are  left  without  next  amplitude,  YOU (2): 

i 

Z„<1)((,)  =  i,x  10'10. 


To  obtain  YOU(3)  through  YOU(Nl)  (or  through  YOU(MATCH)  if  MATCH  <  Nl,  />., 
z„  <  H\)  we  employ  the  "three-point  extrapolation"  finite-difference  equivalent  of  Eq.  (Al).  For 
j  —  2,  ...  ,  Nl  — 1  we  have: 
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zruf,) 


-Z„n,((j  -  2)f,) 


z.u,(0-  D/i) 


24  -  10 


to2f 2 


12  + 


tu2r 2 


ci2  (0-2)/,) 


-  *„2r,2 


c2  (O'  -  1)/,) 

12  + 


-  k2^ 


<v2r,2 


c,2Ofi) 


—  fc„V 


If  MATCH  >  Nl,  the  above  calculation  will  procede  all  the  way  to  the  boundary  value  of  the 
eigenfunction  in  the  water  at  the  water/sediment  interface:  YOU(Nl)  =  Zna,(H i).  From  the  boundary 
conditions,  we  immediately  obtain  the  boundary  value  of  the  eigenfunction  in  the  sediment  layer: 
YOU(NIPLSI)  =  Zn(2)(Nt)  =  (P1/P2)  Z„n)(W,).  To  obtain  the  next  point,  we  again  use  the  Taylor 
polynomial,  this  time  expanded  about  H\. 


Z„,2'(//,  +  t2)  =  Z„(2I(W,)  + 


<2 


dZ ' 


dz 


d2Z 


(2) 

n 


dz 2 


h ' 


Now 


dZ <2’ 

dZ"}  | 

dZ„ 

dz 

dz  1 

dz 

because  of  the  boundary  conditions.  We  program  the  derivative  numerically: 

-  ( — 3Z„<1)(//1  —  2 /,)  +  4 Ziu(H,-  /,)  -  Zi"{Hi))ntx. 
From  Eq.  (Al),  we  have  the  second  derivative: 


dZ„ 

dz 


2  7*2) 


</2Z, 


dz2 


2 

—  _ , 

"  k2 

C22(//>)  1 

Z„t2>  (//,). 


Substituting  these  values  into  the  Taylor  polynomial,  one  obtains  the  value  of  Z„  (//,  +  t2 )• 


From  this  point  on,  i.e.,  for  calculation  of  YOUCN1PLS3)  through  YOU(MATCH),  we  again 
make  use  of  the  finite-difference  equivalent  of  Eq.  (Al);  the  procedure  is  straightforward  and  we  shall 
not  repeat  the  details. 


The  quantity  YOU(MATCH)  is  the  value  of  the  eigenfunction  at  zm,  as  calculated  by  the  original 
procedure  of  integrating  upward  from  the  bottom.  The  variable  PPLUS1  is  the  same  quantity,  but  cal¬ 
culated  by  the  present  procedure  of  starting  at  the  top  and  integrating  downward.  It  remains  to  "match" 
the  solutions,  which  is  accomplished  by  the  DO  149  loop.  Each  of  the  values,  from  the  surface  to  zm, 
is  multiplied  by  the  factor  YOU(MATCH)/PPLUSl.  This  ensures  the  continuity  of  the  eigenfunction 
at  zm.  Continuity  of  its  first  derivative  follows  from  the  condition  that  the  correct  eigenvalue  has 
already  been  determined.  Continuity  of  higher  derivatives  follows  from  Eq.  (Al).  This  completes  the 
task,  since  normalization  takes  place  later  in  SOLID. 

With  control  returned  from  HALFS  to  SOLID,  the  eigenfunction  is  normalized,  as  stated  above. 
Next,  the  scattering  ratios  and  the  bottom  attenuation  ratios  are  determined  (see  Eqs.  (A3)  and  (A4)). 
The  quantity  T0,n  is  termed  AIRH20U)  and  T,  „  is  termed  H2O2ND0),  where  1  designates  the  mode 
order.  The  quantity  y(„2)  is  termed  ATTEN20),  y”r>  is  termed  ATTENC(l),  and  is  termed 
ATTENS(l).  The  values  of  these  parameters  are  as  defined  in  the  old  version  of  the  program  (see  Ref. 
Al);  however,  the  code  has  been  modified  to  take  into  account  the  discus.  '  procedure  of  scaling  the 
eigenfunction.  The  only  numerical  change  is  in  the  formula  used  for  calculating  derivatives  of  the 
eigenfunction  at  the  surface  and  bottom.  For  example,  whereas  we  previously  lAl]  defined 

DYOUA - (YOU(2) »ANORM-YOU(  1 )  *ANORM)/DLl, 
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we  now  use  the  more  accurate  formula: 


DYOUA*=  (— 3.*YOU(l)  *ANORM  +  4.»YOU(2)  »ANORM  -YOU(3)  •ANORM)/2./DLl. 


In  the  present  version  of  SOLID,  normalization  is  accomplished  first,  thus  we  have  programmed 
YOU(J)»ANORM  as  UNRMI(I,J)  or  UNRM2(I,J),  where  the  mode  order  I  has  been  inserted. 

8  =  «2Y*2’  +  +  « Jsr  +  So..  +  Si.,  +  a„.  (A3) 


So  .n 

S\.n 


c,<  0) 

Oi 


Cl 


(//,) 


-  k„ 


O./i 


\.n 


(A4) 


In  many  cases,  attenuation  of  the  modal  field  due  to  absorption  by  the  water  is  not  appreciable. 
In  some  cases,  however,  this  loss  mechanism  has  been  found  to  yield  a  significant  contribution,  thus  it 
is  now  included  in  the  transmission  loss  calculated  by  MOATL  via  Eq.  (A3).  In  a  manner  similar  to 
that  of  the  sediment  layer  and  basement,  we  may  write  the  attenuation  as  a„  =  e^1*,  where  «!  is  the 
plane-wave  absorption  coefficient  in  an  infinite  medium  of  ocean  water.  Whereas  e2,  e3c,  and  e3s 
depend  on  the  particular  sediment  and  basement  chosen  (and  are  input  variables),  e|  may  by  taken  as 
constant  and  given  empirically  by  [A21: 


«i  “ 


(0.1)  f1 

1+/2 


40/2 

4100+/2 


+  2.75  x  1(T4/2, 


where  /  is  the  frequency  in  kHz  and  et  is  given  in  dB/kyd.  The  coefficient  ct  is  termed  EPW  by 
SOLID,  where  it  is  convened  to  nepers/m.  The  quantity  y  ‘n  measures  the  degree  of  interaction  of  the 
nth  mode  with  the  absorption  mechanism  and  is  given  (in  a  form  similar  to  that  of  y,52)),  by: 


yln"  =  Pi 


I  Uff  1  >  (  2  )  ]  2 

Ci(z) 


The  superscript  on  the  eigenfunction  refers  to  the  layer  number,  in  this  case  the  water  layer. 
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Appendix  B 

FORTRAN  NAMELIST  OF  PROGRAM  MOATL  VARIABLES 
AND  SYSTEM  SUBROUTINES  AND  FUNCTIONS 


All  the  variables  of  PROGRAM  MOATL  are  alphabetically  listed  below  with  their  corresponding 
definitions.  For  the  arrays  in  the  list,  the  index  N  denotes  mode-order,  the  index  1  denotes  range,  and 
the  variables  J1  and  J2  denote  receiver  and  source  identification,  respectively.  Occasionally,  reference 
is  made  to  the  "SR  group"  or  the  "LR  group."  For  further  explanation  and  definitions,  see  the  section 
"Step-By-Step  Analysis— Part  3." 

Most  of  the  variables  occurring  in  the  subroutines  have  been  described  elsewhere  [B1  ];  some 
name  changes  and  new  variables  are  defined  in  Appendix  A. 

The  computer  system  library  functions  appearing  in  the  FORTRAN  code  are  described  at  the  end 
of  this  Appendix  in  Table  Bl. 

The  on-line  plotter  software  includes  the  following  subroutines,  which  are  not  discussed  here: 
CENTRE,  CNUMBR,  ENDPLT,  NXAXIS,  NYAXIS,  ORIGIN,  PLOT,  PLOTS,  and  SYMBOL. 

ABSOR(N).  An  array  containing  the  water  absorption,  a„,  calculated  in  the  subroutines. 

AIRH20(N):  An  array  containing  the  scattering  ratios  r0.„,  used  in  calculating  the  attenuation  due  to 
interaction  of  the  various  modes  with  a  statistically  rough  boundary  at  the  air/water 
interface. 

AKN(N):  An  array  containing  the  eigenvalues  (modal  wave  numbers),  k„ ,  as  calculated  by  the  subrou¬ 
tines. 

ATTENC(N):  An  array  containing  the  basement  compressional  attenuation  ratios,  y^3l). 

ATTENS(N):  An  array  containing  the  basement  shear  attenuation  ratios,  yli5). 

ATTEN2(N):  An  array  containing  the  sediment  attenuation  ratios,  y(n2). 

Al:  The  number  plus  fraction  (starting  at  the  air/water  surface)  of  the  incremental  layer  corresponding 
to  the  source  or  receiver  depth. 

A2:  In  a  range-dependent  model,  Al  is  used  (as  defined  above)  for  the  SR  calculation  and  A2  (defined 
the  same  way)  is  used  for  the  LR  calculation. 

CBS:  the  range-interpolated  value  of  c, (//) ,  the  sound  speed  at  the  bottom  of  the  water  layer. 

CB1:  The  value  of  C\(H)  at  the  SR  range. 
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CB2:  The  value  of  ct(H)  at  the  LR  range. 

COMP:  The  compressional  velocity,  c 3(.,  of  the  basement. 

COPL(Jl,J2,l):  An  array  containing  the  coherent  transmission  loss  or  transmission-loss  anomaly  in 
decibels. 

CTS:  The  range-interpolated  value  of  c,(0),  the  sound  speed  at  the  surface  of  the  water  layer. 

CT1:  The  value  of  ci (0)  at  the  SR  range. 

CT2:  The  value  of  c^O)  at  the  LR  range. 

C 1  ( J ) :  An  array  containing  the  sound-speed  profile  in  the  water  layer,  ci(z). 

C2(J):  An  array  containing  the  sound-speed  profile  in  the  sediment  layer,  c2(z). 

DBMAX:  The  maximum  loss  (in  dB)  on  the  vertical  axis  of  the  plot. 

DBMIN:  The  minimum  loss  (in  dB)  on  the  vertical  axis  of  the  plot. 

DEPRE(Jl):  An  array  containing  the  receiver  depths. 

DEPSOU2):  An  array  containing  the  source  depths. 

DLTA1:  That  fraction  of  an  incremental  layer  by  which  the  source  or  receiver  depth  exceeds  the  depth 
of  the  nearest  shallower  incremental  layer  boundary. 

DLTA2:  In  a  range-dependent  model,  DLTA1  is  used  (as  defined  above)  for  the  SR  calculation  and 
DLTA2  (defined  the  same  way)  is  used  for  the  LR  calculation. 

DR:  The  distance  (increment)  between  ranges  at  which  transmission  loss  (or  loss  anomaly)  is  calcu¬ 
lated. 

DX:  The  number  of  kilometers  per  inch  on  the  horizontal  (range)  axis  of  the  plot. 

DY.  The  number  of  decibels  per  inch  on  the  vertical  (loss)  axis  of  the  plot. 

EIGVLl(N):  An  array  containing  the  eigenvalues,  k„ ,  at  the  SR  range. 

EIGVL2(N):  An  array  containing  the  eigenvalues,  k„ ,  at  the  LR  range. 

EPSLN.  The  criterion  for  the  accuracy  of  the  determined  eigenvalues  and  eigenfunctions,  i.e.,  the 
amount  by  which  the  air/ water  surface  value  of  the  normalized  eigenfunction  may 
differ  from  the  pressure-release  boundary  condition  requirement  of  being  identically 
zero. 

EP1:  The  plane-wave  absorption  coefficient  of  the  sediment,  e2,  expressed  in  dB/Hz-m. 

EP2:  The  basement  compressional  plane-wave  absorption  coefficient,  e3c,  expressed  in 
dB/Hz-m. 

EP3:  The  basement  shear  plane-wave  absorption  coefficient,  e3(,  expressed  in  dB/Hz-m. 
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EP4:  Same  as  EP1,  expressed  in  nepers/m. 

EP5:  Same  as  EP2,  expressed  in  nepers/m. 

EP6:  Same  as  EP3,  expressed  in  nepers/m. 

F:  The  source  frequency,  f. 

GA(N):  An  array  containing  the  range-interpolated  values  of  the  water  absorption,  a„. 

GB(N):  An  array  containing  the  range-interpolated  values  of  the  scattering  ratios  T  | 

GT(N):  An  array  containing  the  range-interpolated  values  of  the  scattering  ratios  r0„. 

GUN):  An  array  containing  the  range-interpolated  values  of  the  attenuation  ratios  y{„2\ 

G2(N):  An  array  containing  the  range-interpolated  values  of  the  attenuation  ratios 

G3(N):  An  array  containing  the  range-interpolated  values  of  the  attenuation  ratios  ■y^s’. 

HS.  The  range-interpolated  value  of  the  water  layer  thickness  (»>.,  the  water  depth),  //,. 

H10:  The  value  of  the  water  layer  thickness  H\  at  the  receiver,  i.e.,  at  zero  range. 

HI  1:  The  value  of  the  water  layer  thickness  H\  at  the  SR  range. 

HI 2:  The  value  of  the  water  layer  thickness  at  the  LR  range. 

H2:  The  value  of  the  sediment  layer  thickness  H2- 

H202ND(N):  An  array  containing  the  scattering  ratios  Fj.,,,  used  in  calculating  the  attenuation  due  to 
interaction  of  the  various  modes  with  a  statistically  rough  boundary  at  the 
water/sediment  interface. 

1.  A  DO-LOOP  index  used  primarily  for  the  DO  300  range  loop. 

I A 1 :  The  truncated  integer  value  of  Al. 

IA2:  The  truncated  integer  value  of  A2. 

1FREQ:  The  DO  340  frequency  loop  index. 

Ill:  An  input  flag  used  to  specify  whether  the  basement  is  to  be  modeled  as  a  fluid  or  a  solid. 

IM:  A  DO-LOOP  index  denoting  mode-order. 

INC1:  A  variable  having  a  value  one  greater  than  the  number  of  depth-increment  values  to  be  skipped 
between  printed-out  values  of  the  eigenfunctions  in  the  water  layer. 

INC2:  A  variable  having  a  value  one  greater  than  the  number  of  depth-increment  values  to  be  skipped 
between  printed-out  values  of  the  eigenfunctions  in  the  sediment  layer. 

IPLOT:  An  input  flag  used  to  specify  whether  or  not  plotting  of  transmission  loss  (or  loss  anomaly)  vs 
range  is  to  be  included  in  the  output. 
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ISPRD:  An  input  flag  used  to  specify  whether  transmission  toss  or  transmission  loss  anomaly  is  to  be 
calculated. 

J.  An  implied  DO-LOOP  index  used  for  the  sound-speed  profiles. 

J 1 :  A  DO-LOOP  index  used  for  denoting  the  various  receiver  depths. 

J2:  A  DO-LOOP  index  used  for  denoting  the  various  source  depths. 

KA:  A  parameter  whose  value  flags  the  normal-mode  subroutines  to  store  the  calculated  eigenfunc¬ 
tions  in  either  the  SR  array  UNRM1(N,K)  or  the  LR  array  UNRM2(N,K). 

LI1:  The  number  of  incremental  steps  into  which  the  water  layer  is  to  be  divided  for  the  numerical  cal¬ 
culations. 

L12:  The  number  of  incremental  steps  into  which  the  sediment  layer  is  to  be  divided  for  the  numerical 
calculations. 

MDPRNT:  An  input  flag  used  to  specify  whether  or  not  calculated  normal-mode  amolitude  functions 
(eigenfunctions),  u„(z),  are  to  be  printed  out. 

MGNTD:  A  machine-dependent  parameter  set  equal  to  approximately  half  (but  not  larger  than  half) 
the  dynamic  range  of  a  single-precision  real  constant. 

NDRE:  The  number  of  receiver  depths  for  which  calculations  are  to  be  performed. 

NDSO:  The  number  of  source  depths  for  which  calculations  are  to  be  performed. 

ND1:  The  number  of  points  (depths)  in  the  water  layer  at  which  the  sound  speed  is  supplied. 

ND2:  The  number  of  points  (depths)  in  the  sediment  layer  at  which  the  sound  speed  is  supplied. 

NMFREQ:  The  number  of  source  frequencies  (or  more  accurately,  the  number  of  separate  environ¬ 
mental  cases)  for  which  the  entire  transmission-loss  calculations  are  to  be  performed. 

NMODE.  The  total  number  of  modes  that  exist  and/or  that  are  to  be  used  in  the  calculations. 

NRBUF:  The  number  of  range  points  at  which  environmental  data  are  to  be  supplied. 

NRCALC:  The  number  of  range  points  at  which  loss  calculations  are  to  be  performed. 

Nl:  LI  1  +  1  in  the  subroutines. 

Nil:  LI1  + 1  at  the  SR  environmental  data  set. 

N12:  HI  +  1  at  the  LR  environmental  data  set. 

N2:  LI2+1  in  the  subroutines. 

PHASE(J1,J2):  An  array  containing  the  phase  angles  in  degrees  of  the  coherent  transmission  loss. 

PL  (J  1 ,  J2,I) :  An  array  containing  tne  incoherent  transmission  loss  or  transmission-loss  anomaly  in  deci¬ 
bels. 
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PLTAR(K):  An  array  used  solely  by  the  machine  for  buffering  the  execution-generated  plot  commands 
onto  Calcomp  plotter  tape. 

PMGTD:  The  value  of  ten  raised  to  the  power  MGNTD. 


PPRCN.  The  value  of  ten  raised  to  the  power  (—  PRCSN). 

PRCSN:  A  machine-dependent  parameter  specifying  one  less  than  the  approximate  number  of  decimal 
digits  associated  with  the  mantissa  of  a  DOUBLE  PRECISION  real  constant. 


Q.  A  variable  containing  the  product  of  Q1  and  Q2. 


QC(J1,J2): 


An  array  containing  the  intermediate  results,  £  — 

n=- 1 

1  ^  AT  <  A,  which  are  used  in  calculating  the  coherent  loss. 


(lp)u„(l)rh 

i  [h 


-  A  r 

e 


cos  <(>„ , 


QQ:  The  average  attenuation  coefficient,  A„,  times  the  range,  r. 


QS  ( J 1 ,  J  2) .  An  array  containing  the  intermediate  results,  £ 


,  u„U0)u'„(Or  -a 


n-i 


,  V: 

<Pn 


e  "  sin  <t> 


Ql:  The  term 


I  <  A'  <  A,  which  are  used  in  calculating  the  coherent  loss. 
uMo)u'„(Qr'h 

<t>n 


Q2:  The  term  e  " . 

RAN(I):  An  array  containing  the  range  points,  r,  at  which  results  are  to  be  calculated. 


RANGEl:  The  range  of  the  SR  environmental  data  set. 

RANGE2:  The  range  of  the  LR  environmental  data  set. 

RA1  (N):  An  array  containing  the  water  absorption,  a„  at  the  SR  range. 

RA2(N):  An  array  containing  the  water  absorption,  a„,  at  the  LR  range. 

RBl(N):  An  array  containing  the  scattering  ratios  T  j  „  at  the  SR  range. 

RB2(N):  An  array  containing  the  scattering  ratios  Ti.,  at  the  LR  range. 

RE(J1,N):  An  array  containing  the  depth-interpolated  values  of  w„(£0 )«  which  are  the  eigenfunction 
values  at  the  various  receiver  depths. 

REC:  A  PARAMETER  which  specifies  the  maximum  number  of  receiver  depths. 

REC5:  A  PARAMETER  which  is  greater  than  or  equal  to  REC  and  which  is  an  integral  multiple  of 
five. 

RERHOl:  The  water-layer  density  p,  at  the  receiver,  i.e.,  at  zero  range. 

RERH02:  The  ratio  of  the  density  of  the  sediment  layer  to  that  of  water,  i.e.,  p2,  at  the  receiver  range. 
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RHOl:  The  water-layer  density  p^  at  the  variable  range  r. 

RH02:  The  ratio  of  the  density  of  the  sediment  layer  to  that  of  water,  i.e.,  p2,  at  the  variable  range  r. 

RH03:  The  ratio  of  the  density  of  the  basement  to  that  of  water,  i.e.,  p 3,  at  the  variable  range  r. 

RMAX:  The  range  of  the  final  input  environmental  data  set  and/or  the  maximum  range  at  which  cal¬ 
culations  are  desired. 

RMXKM:  The  length,  in  km,  of  the  range  axis  for  plotted  output. 

RNG.  A  PARAMETER  which  specifies  the  maximum  number  of  range  points  at  which  loss  calcula¬ 
tions  are  to  be  performed. 

RT1(N):  An  array  containing  the  scattering  ratios  rn  „  at  the  SR  range. 

RT2(N):  An  array  containing  the  scattering  ratios  r0  „  at  the  LR  range. 

RU(N):  An  array  containing  the  sediment-attenuation  ratios  y{n2)  at  the  SR  range. 

R12(N):  An  array  containing  the  sediment-attenuation  ratios  yi2}  at  the  LR  range. 

R21(N):  An  array  containing  the  basement  compressiona!  attenuation  ratios  y  ',3‘ '  at  the  SR  range. 

R22(N):  An  array  containing  the  basement  compressional  attenuation  ratios  -y,'3‘ 1  at  the  LR  range. 

R 3 1  ( N ) :  An  arrav  containing  the  basement  shear  attenuation  ratios  y,1,3'1  aI  the  SR  range. 

R32(N):  An  array  containing  the  basement  shear  attenuation  ratios  y!,u>  at  the  LR  range. 

SA(N)  An  array  containing  the  values  of  the  integrals  f  a„(r )  dr. 

SB(N):  An  array  containing'the  values  of  the  integrals  f  T.  „(r)  dr. 

SCALE:  A  variable  which  contains  the  interpolation  factor  A  used  in  the  linear  range  interpolation. 

(See  the  section  "Step-By-Step  Analysis— Part  7,”  specifically  Eqs.  (16)  and  (17).) 

SE(N):  An  array  containing  the  values  of  the  integrals  </>„  *  J*n  fc„(r)  dr. 

SHEAR:  The  shear  velocity,  c3s,  of  the  basement. 

S1G0:  The  root-mean-square  wave  height  o-0,  used  in  calculating  the  term  S0  n. 

SIG1:  The  root-mean-square  excursion  of  the  water/sediment  interface  tr,  ,  used  in  calculating  the 
term  A,  „. 

SM(J2,N):  An  array  containing  the  depth-interpolated  values  of  w'„(£),  which  are  the  eigenfunction 
values  at  the  various  source  depths. 

SOC:  A  PARAMETER  which  specifies  the  maximum  number  of  source  depths. 

ST(N):  An  array  containing  the  values  of  the  integrals  ;;  r„.„(r)  dr. 
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S1(N):  An  array  containing  the  values  of  the  integrals  f  y l2)(r)  dr. 

S2(N):  An  array  containing  the  values  of  the  integrals  J"o  y,!3‘  '(/■)  dr. 

S3(N):  An  array  containing  the  values  of  the  integrals  f  y  {„is)(r)  dr. 

%/o 

TH:  The  name  of  a  COMMON  block  used  for  variables  common  to  PROGRAM  MOATL,  SUBROU¬ 
TINE  HALFF,  and  SUBROUTINE  HALFS. 

TITLE(K):  Any  alphanumeric  label  used  to  identify  the  computer  run. 

TN:  The  name  of  a  COMMON  block  used  for  variables  common  to  PROGRAM  MOATL,  SUBROU¬ 
TINE  FLUID,  and  SUBROUTINE  SOLID. 

TNH:  The  name  of  a  COMMON  block  used  for  variables  common  to  PROGRAM  MOATL,  SUBROU¬ 
TINE  FLUID,  SUBROUTINE  SOLID,  SUBROUTINE  HALFF,  and  SUBROUTINE 
HALF  S. 

TNI:  The  name  of  a  COMMON  block  used  for  '’ariables  common  to  PROGRAM  MOATL,  SUBROU¬ 
TINE  FLUID,  SUBROUTINE  SOLID.  SUBROUTINE  ITRTF,  and  SUBROUTINE 
ITRTS. 

TNIH.  The  name  of  a  COMMON  block  used  for  variables  common  to  PROGRAM  MOATL  and  all 
subroutines. 

UNRM1  (N,J2):  An  array  containing  the  eigenfunctions  u,j(z)  at  the  SR  range. 

UNRM2(N,J2):  An  array  containing  the  eigenfunctions  u'n(z)  at  the  LR  range. 

WN(N):  An  array  containing  the  range-interpolated  values  of  the  eigenvalues  k„. 

X:  The  x-coordinate  of  a  point  to  be  plotted  on  the  graph. 

XLENG:  The  to-scale  length,  in  inches,  of  the  range  axis  for  the  plotted  output  graph. 

XXRl.  The  depth-interpolated  value  of  the  eigenfunction  at  the  source  depth,  at  the  SR  range. 

XXR2:  The  depth-interpolated  value  of  the  eigenfunction  at  the  source  depth,  at  the  LR  range. 

X2:  XLENG  divided  by  two. 

Y:  The  y-coordinate  of  a  point  to  be  plotted  on  the  graph. 

YLENG:  The  length,  in  inches,  of  the  loss  axis  for  the  plotted  output  graph. 

Y2:  YLENG  divided  by  two. 

Z1  (J):  An  array  containing  the  sound-speed-profile  depth  points  in  the  water  layer. 

Z2(J):  An  array  containing  the  sound-speed-profile  depth  points  in  the  sediment  layer. 
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Table  B1  —  System  Library  Functions 


Form 

Definition 

Mode 

of 

Argument 

Mode 

of 

Result 

ABS(X) 

Absolute  value  of  X 

Real 

Real 

ALOGIO(X) 

Logarithm  to  base  10  of  X 

Real 

Real 

ATAN2(X,Y) 

Arctangentt  of  (X  divided  by  Y) 

Real 

Real 

DCOS(D) 

Cosine  of  D+ 

Double 

Double 

DSIN(D) 

Sine  of  D+ 

Double 

Double 

DSQRT(D) 

Square  root  of  D 

Double 

Double 

EXP(X) 

e  raised  to  power  X 

Real 

Real 

FLOAT(I) 

Convert  integer  1  to  real 

Integer 

Real 

1F1X(X) 

Truncate  real  X  to  integer 

Real 

Integer 

tin  radians 


REFERENCES 

BI.  J.F.  Miller  and  F.  Ingenito,  NRL  Memorandum  Report  3071,  June  1975. 
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FORTRAN  LISTING  AND  CROSS-REFERENCE 
OF  PROGRAM  AND  SUBROUTINES 

The  FORTRAN  code  for  each  of  the  main  program  and  subroutines  is  listed  on  the  following 
pages.  Following  each  listing  is  a  cross-reference  for  that  routine,  which  gives  all  of  the  variable  and 
routine  names  used,  in  alphabetical  order.  Following  each  FORTRAN  name  in  the  cross-reference  is  a 
list  of  control  statement  numbers  (CS Ns)  referencing  the  statements  in  which  the  name  appears. 

The  following  control  statements  are  indented  three  spaces  for  ease  in  following  the  FORTRAN 
code:  (1)  DO  loop,  (2)  GO  TO  statement,  (3)  transfer-of-control  IF  statement,  and  (4)  calls  to  subrou¬ 
tines. 

Each  Part  of  PROGRAM  MOATL  (as  subdivided  for  the  discussion  in  the  "Step-By-Step 
Analysis”  section)  is  marked  with  a  COMMENT. 
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